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ABSTRACT 

We report the results of an investigation of particle acceleration and electron-positron plasma 
generation at low altitude in the polar magnetic flux tubes of rotation-powered pulsars, when 
the stellar surface is free to emit whatever charges and currents are demanded by the force- 
free magnetosphere. We apply a new ID hybrid plasma simulation code to the dynamical 
problem, using Particle-in-Cell methods for the dynamics of the charged particles, including a 
determination of the collective electrostatic fluctuations in the plasma, combined with a Monte 
Carlo treatment of the high-energy gamma-rays that mediate the formation of the electron- 
positron pairs. We assume the electric current flowing through the pair creation zone is fixed by 
the much higher inductance magnetosphere, and adopt the results of force-free magnetosphere 
models to provide the currents which must be carried by the accelerator. The models are 
spatially one dimensional, and designed to explore the physics, although of practical relevance 
to young, high-voltage pulsars. 

We observe novel behaviour (a) When the current density j is less than the Goldreich-Julian 
value (0 < j/joj < 1), space charge limited acceleration of the current carrying beam is mild, 
with the full Goldreich-Julian charge density comprising the charge densities of the beam and 
a cloud of electrically trapped particles with the same sign of charge as the beam. The voltage 
drops are of the order of mc 2 /e, and pair creation is absent, (b) When the current density 
exceeds the Goldreich-Julian value (//j’gj > 1), the system develops high voltage drops (TV 
or greater), causing emission of curvature gamma-rays and intense bursts of pair creation. The 
bursts exhibit limit cycle behaviour, with characteristic time-scales somewhat longer than the 
relativistic fly-by time over distances comparable to the polar cap diameter (microseconds), (c) 
In return current regions, where j/jos < 0, the system develops similar bursts of pair creation. 
These discharges are similar to those encountered in previous calculations by Timokhin of pair 
creation when the surface has a high work function and cannot freely emit charge. In cases 
(b) and (c), the intermittently generated pairs allow the system to simultaneously carry the 
magnetospherically prescribed currents and adjust the charge density and average electric field 
to force-free conditions. We also elucidate the conditions for pair creating beam flow to be 
steady (stationary with small fluctuations in the rotating frame), finding that such steady flows 
can occupy only a small fraction of the current density parameter space exhibited by the force- 
free magnetospheric model. The generic polar flow dynamics and pair creation are strongly 
time dependent. The model has an essential difference from almost all previous quantitative 
studies, in that we sought the accelerating voltage (with pair creation, when the voltage drops 
are sufficiently large; without, when they are small) as a function of the applied current. 
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The ID results described here characterize the dependence of acceleration and pair creation 
on the magnitude and sign of current. The dependence on the spatial distribution of the current 
is a multi-dimensional problem, possibly exhibiting more chaotic behaviour. We briefly outline 
possible relations of the electric field fluctuations observed in the polar flows (both with and 
without pair creation discharges) to direct emission of radio waves, as well as revive the 
possible relation of the observed limit cycle behaviour to microstructure in the radio emission. 
Actually modelling these effects requires the multi-dimensional treatment, to be reported in a 
later paper. 

Key words: acceleration of particles -plasmas -stars: magnetic field -stars: neutron - 
pulsars: general. 


1 INTRODUCTION 

Young pulsar wind nebulae (PWNe) show that rotation-powered 
pulsars (RPPs) have dense magnetospheres, at least with regard 
to those regions that feed the plasma outflow (e.g. Bucciantini, 
Arons & Amato 2011). Electron-positron pair creation in the open 
field line region that connects to the external world is the only 
known candidate for the origin of such outflows, with acceleration 
and convertible gamma-ray emission occurring either at low alti- 
tude (Sturrock 1971) or in the outer magnetosphere (Cheng, Ho & 
Ruderman 1986). High-density flows that can feed all the open field 
lines can exist only in the low-altitude polar cap region (for a general 
discussion, see Arons 2009). 

Theoretical studies of charged particle flow from the magnetic 
polar regions of RPPs began with the observation by Goldreich 
& Julian (1969) that an isolated magnetic rotator in vacuum must 
have a charged magnetosphere almost corotating with the star. Since 
RPPs are strongly gravitationally bound and cool objects (thermal 
scale height in any atmosphere orders of magnitude less than the 
stellar radius) and have no external source of plasma supply (so far 
as we know), the only plasma source is an extraction of charged 
particles from the stellar surface, leading to a conjectured magneto- 
sphere whose plasma is fully charge separated, in contrast to all other 
known astrophysical systems, whose plasmas are charged but quasi- 
neutral. Goldreich & Julian (1969) speculated that on polar field 
lines - those that extend beyond the light cylinder located at cylin- 
drical radius R L c = cP/2n — 48 000 P km, P = rotation period in 
seconds - a charge-separated outflow would form. They argued that 
the energy/particle in the outflow would be no more than the gravi- 
tational escape energy GM t /R t ~ 0.3mc 2 (M*/1.4MQ)(10km/R*); 
M* and R* are star’s mass and radius correspondingly - the par- 
ticles leave at non-relativistic energies in spite of the fact that 
the electric potential drop across the polar field lines is equal 
to the full potential of an open rotating magnetosphere with a 
dipole magnetic field <t> m = ■dW^Jc ~ lOfP/lO -15 ) 1 ^ 2 ^^ 3 ^ 2 TV, 
W R = rotational energy loss rate, with P = dP/df. <f> m vastly 
exceeds the rest energy and gravitational energy of the particles, 
either electrons or protons (or He or other ions populating the 
star’s crust and atmosphere). The super strong magnetic field sup- 
presses free acceleration of the particles in the transverse electric 
field, whose primary (‘zeroth order’) consequence is the corota- 
tion of the field lines with the magnetic field embedded in the 
neutron star (NS), with the field line motion measured by the 
E x B drift of charged particles across the magnetic field (which 
occurs even when the particles have zero Larmor gyration). The 
particle loss rate in the conjectured charge-separated scenario is 
N R = c<t> m /e i» 2 x 10 30 (P/10- 15 ) 1 / 2 P- 3 '' 2 s -1 , orders of magni- 
tude less than that inferred from the injection of plasma into the 
young PWNe. The electrodynamics of the magnetosphere differ 


drastically, depending on whether the particle loss rate falls short 
of or exceeds N R . For the young PWNe, the particle injection rate 
exceeds N R (by a lot). In that case, the magnetosphere’s basic state 
should be one in which E ■ B = 0, with no parallel acceleration 
sufficient to generate convertible gamma-rays occurring under the 
pulsar’s rotational control. 

The discovery of gamma-ray pulsars in the 1970s, and their pro- 
liferation into a population with more than 100 such stars in the 
most recent published Fermi pulsar catalogue (Abdo et al. 2010), 
has shown that parallel acceleration to GeV gamma-ray emitting 
energies (indeed, multi-hundred GeV, in the Crab pulsar; Aliu et al. 
2011) must occur somewhere, with energy efficiency exceeding a 
few tenths of a per cent, as measured by i y /W R , L y = gamma-ray 
luminosity. If the acceleration is limited by radiation reaction, as 
is true in many models, L y is a good proxy for the energy put into 
parallel acceleration. L y /W R can approach as much as 50 per cent 
at smaller spin down luminosities. Just how some fraction of the to- 
tal potential drop gets released in acceleration along B, gamma-ray 
emission and pair creation has been mysterious since the beginning 
of pulsar research, made relevant to the real world by the gamma-ray 
discoveries. Since the polar cap source is the only one capable of 
feeding the whole (open) magnetosphere, its understanding remains 
of central interest to modelling pulsar magnetospheres, even though 
the spectral and beaming characteristics of the pulsed gamma-rays 
are better modelled by accelerators in the outer magnetosphere. 

Free particle outflow from the NS surface is a common assump- 
tion in most of the current pulsar models (see Section 2). The polar 
cap accelerator problem has been studied under that assumption be- 
fore (e.g. Michel 1974; Fawley, Arons & Scharlemann 1977; Mestel 
et al. 1985; Shibata 1997; Beloborodov 2008). Michel (1974) and 
Fawley et al. (1977) obtained solutions for the non-neutral space 
charge limited charged particle flow for the current density almost 
equal to the GJ current density. All these models assume strictly 
steady flow in the corotating frame on all time-scales. In these 
models, the charge density of the current carrying beam supplies 
almost all of the charge density needed to short out the parallel 
component of the electric field, while leaving a residuum E\ \ suffi- 
cient to accelerate the beam - relativistic energies in a temporally 
steady flow are found if the current density y)| = — (B/P) cos x + 
small corrections = jc j; X — ^(P- H), the pulsar inclination angle. 
Mestel et al. (1985) showed that the velocity of the beam is mono- 
tonically increasing with altitude to relativistic speeds only if the 
current density is larger than Jqj . If the current density is smaller than 
ja the temporally steady velocity of the beam (assumed to have no 
momentum dispersion) oscillates spatially, i.e. particles accelerate 
and decelerate to a complete halt as they move outwards into the 
magnetosphere. Beloborodov (2008) rediscovered Mestel et al.’s 
solution and suggested that in the region of the polar cap where 
III < Jgj particles will not be accelerated up to high energies as 
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the beam velocity oscillates, but along the magnetic field lines with 
7 II / /gj > 0 or 7 1 1 /7 gj < 0 particle acceleration will be efficient and will 
lead to pair formation. The quantitative model which we describe 
in this paper lends some support to Beloborodov’s speculations, 
although it does not agree with them in detail. 

In this paper we describe our study of the physics of the polar cap 
accelerator in the space charge limited flow regime starting from 
first principles - assuming free particle outflow from the surface of a 
NS we compute the electric field, particle acceleration, gamma-ray 
emission, propagation and pair creation simultaneously. It extends 
the study of current flow and pair cascades in NS magnetospheres 
using the theoretical formulation and self-consistent numerical tech- 
niques introduced in Timokhin (2010). 

The plan of the paper is as follows. In Section 2, we review the 
properties of the current flow imposed by the magnetosphere, in the 
force-free model, pointing out that the current density is the main 
parameter which regulates the efficiency of particle acceleration. In 
Section 3 we review the properties of stationary solutions for the 
charge-separated space charge limited flow problem. In Section 4 
we briefly describe our numerical model. In Section 5 we describe 
the results of numerical modelling for the case of sub-GJ current 
density, the regime when particle acceleration is inefficient and no 
pair creation is possible. In Section 6 we consider flow regimes with 
efficient pair creation. In Section 6.4 we pay special attention to the 
stationary flow regime, which up to now was assumed in (most) 
works on pulsar polar cap accelerators, and describe why it has lim- 
ited relevance to the force-free model of the pulsar magnetosphere. 
We discuss the implications of our results for the physics of RPPs 
in Section 7 and summarize our conclusions in Section 8. 

2 CURRENT DENSITY IN THE POLAR CAP 

Nebular observations of plasma supply by RPPs suggest the open 
field line regions are ‘magnetohydrodynamic (MHD)-like’, i.e. hav- 
ing £|| = 0 except in special zones (such as the polar cap), which 
are, in effect, boundary layers. There is essentially no observational 
information on the properties of the closed magnetosphere, thus the 
simplest hypothesis is to follow Goldreich & Julian (1969) and as- 
sume the magnetosphere is (almost everywhere) filled with plasma 
that shorts out E\ . The magnetosphere has open and closed magnetic 
field line zones. In the open zone plasma flows away into the pulsar 
wind; currents and their associated electromagnetic inertia keep the 
magnetic field open, so sustaining the flow. The plasma flows all 
the way from the base of the open field line zone in the polar cap of 
the pulsar where it is either extracted from the NS surface or gener- 
ated by electron-positron cascades (or both). The distribution of the 
current density across the open field line zone, and therefore across 
the polar cap, is determined by the global magnetospheric structure. 
Stability of the pulsar mean profiles and sharpness of the peaks in 
the spectra of gamma-ray pulsars strongly suggest that on scales 
comparable to the light cylinder, the magnetosphere is stationary 
in the frame corotating with the NS with smooth and continuous 
plasma outflow. However, the stationary corotating magnetosphere 
hypothesis demands stationarity only in a statistical sense; fast lo- 
cal fluctuations which average to a stationary state (and even more 
broadly, global variations) can be included within this picture, so 
long as they do not smear the beaming profiles. 

The polar cap acceleration and possible pair cascade zone - the 
main place where electron-positron plasma that feeds the wind can 
be produced - are much smaller than the characteristic scale R L c of 
the magnetosphere. Therefore its inductance is negligible compared 
with that of the magnetosphere, and the polar current flow within 


the polar cap region must have an average equal to that set by the 
magnetosphere’s global structure. 

In this paper we solve a local problem of how the polar cap 
cascade zone adjusts itself to the current density required by the 
magnetosphere. On the dynamical time-scales typical for the cas- 
cade zone (microseconds) the magnetospherically required current 
density is stationary because it could change only on magneto- 
spheric time-scales (tens of milliseconds up to several seconds). 
The idea that the acceleration zone has a magnetospherically deter- 
mined current appears in the electrodynamics through the magnetic 
induction equation 

3 E« 

— = -47t7j| + c(V x B)|, « -4n(j {] - y m ), (1) 

ot 

with 

c 

jm = ( Y X li magnetosphere )|| (2) 

being the current that sustains the twist to the field lines. In this 
paper we neglect the induced variations in the magnetic field that 
accompany variable £j| - because of the very strong background 
magnetic field (which in the region of interest has V x B = 0), 
these have little dynamical effect on the acceleration and cascade 
dynamics [see Appendix A for the derivation of equation (1)]. 

We study the behaviour of the cascade zone under different cur- 
rent loads, sampling the range of possibilities illustrated in Fig. 1. 
The model has an essential difference from previous studies, in that 
we seek the accelerating voltage (with pair creation when the volt- 
age drops are sufficiently large, without when they are small) as a 
function of the applied current. Previous work has almost entirely 
focused on the opposite direction, seeking the current that emerges 
from the accelerator when the voltage is fixed, either by the ge- 
ometry or by the poisoning of the accelerator by pair creation. In 
addition, we allow the system to be fully time dependent. These 
generalizations lead to qualitatively different results from what has 
appeared before. Our model is one-dimensional, with spatial axis 
along magnetic field lines; from here on we drop subscript || from 
all quantities. 

The characteristic charge density, the Goldreich-Julian charge 
density. 


sets the characteristic current density 

jas = m c ■ (4) 

Following Fawley et al. (1977) and Arons & Scharlemann (1979), 
we also assume that if the star lacks an atmosphere, the work func- 
tion for charged particles to leave the NS surface is small enough 
that any number of charged particles can be extracted from the NS 
surface until the extracting electric field is screened. In the classical 
pulsar regime, the theory of charged particle binding to the crust 
suggests such free emission to be likely (Medin & Lai 2010). More 
relevantly, X-ray observations of heated polar caps (Bogdanov, 
Rybicki & Grindlay 2007, and references therein) suggest that these 
stars have atmospheres overlying the solid and ocean components 
of the crust, which guarantees free emission of charges. The op- 
posite case, with complete suppression of particle emission from 
the surface, was studied in Timokhin (2010), a realization of the 
scenario conjectured by Ruderman & Sutherland (1975). 

As we will show in the following sections there are three qualita- 
tively different plasma flow types in the polar cap of pulsar depend- 
ing on the ratio of the current density imposed by the magnetosphere 
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Figure 1 . Field-aligned current density at the polar cap of the force-free rotator, with / measured in units of the Goldreich— Julian current density /'gj = 
—SI ■ B/Inc. The black circle is the rim of the polar cap - the footprints of the field lines that pass outside the light cylinder fall with that circle. The 
distributed current is shown. The current sheet component coincides with the polar cap boundary. This plot was made by Xuening Bai using results of force-free 
magnetosphere simulations presented in Bai & Spitkovsky (2010). 


to the GJ current density: (i) j m has the same sign and its absolute 
value is smallerthan the GJ current density, 0 < / m //GJ < 1. hereafter 
sub-GJ current density; (ii) j m has the same sign and its absolute 
value is larger than the GJ current density, jm/joj > 1. hereafter 
super-GJ current density; (iii) j m has the opposite sign to the GJ 
current density, jm/jw < 0, hereafter anti-GJ current density. 

The advent of quantitative solutions for the structure of the force- 
free model (e.g. Contopoulos, Kazanas & Fendt 1999; Timokhin 
2006; Spitkovsky 2006; Kalapotharakos & Contopoulos 2009; Bai 
& Spitkovsky 2010) has provided, for the first time, a theory of 
the current flow expected as a function of pulsar parameters ( B , 
P, /, when the magnetic field is a star-centred dipole). Earlier 
modelling of polar cap, slot gap and outer gap accelerators adopted 
the expectation that the current density is of the order of ja, and 
expressed the hope that the accelerator and pair creation physics do 
not sensitively depend on the precise value and spatial distribution 
of j. The results reported here show that the magnitude and sign of 
the current flow do lead to drastic differences in the open field line 
accelerator’s behaviour in the three regimes (i)— (iii), even though 
the order of magnitude of the current is as expected. Fig. 1 , showing 
J = j / ja j, reveals that all three flow regimes occur in the force- 
free magnetosphere model. While | J | always has numerical values 
of the order of unity, it can be negative (return current) as well as 
lying in the separate regimes 0 < j < 1 and J > 1 . We show these 
separate regimes have different dynamical behaviour and different 
implications for pair creation. 

We also show that once the constraints of the steady flow models 
for space charge limited flow are relaxed, the small departures of 
the GJ charge density from the simple estimate —S2B cosx/27Tc 
created by geometric and general relativistic considerations (Arons 
& Scharlemann 1979; Muslimov & Tsygan 1992) that play an es- 
sential role in the steady flow models 1 have little significance when 
the flow is fully time dependent. Since we consider only the polar 
cap region, with altitudes not exceeding the width of the polar flux 
tube r pc = R ts / R. t / R L q = 0.145P -1 ^ 2 km <$( R t = 10 km, spatial 
variation of B is mostly unimportant. If we do not say so explic- 
itly otherwise, throughout the paper we assume that the GJ charge 
density is constant, independent of the distance along B. 

1 If the GJ charge density were uniform and the beam is everywhere rel- 

ativistic and time stationary, the unique model is no acceleration at all 
(Tademaru 1973; Fawley et al. 1977), a severe contradiction. 


3 TIME-STATIONARY SPACE CHARGE 
LIMITED FLOW 

Copious pair creation occurs when there is a sufficiently large ac- 
celerating potential difference along B when pairs are absent. Such 
potential drops, typically > 10 12 V (Sturrock 1971), readily exist 
in the absence of current flow, that is, in a vacuum. Under pulsar 
conditions, if current does flow, pairs appear when the current flow 
co-exists with TV potential drops, 2 that is the current is a relativistic 
beam. Therefore, the starting place is the properties of time station- 
ary, space charge limited, charge-separated flow. In this section we 
give an overview of these properties while detailed derivation of 
the equations used in this section is given in Appendix B together 
with some useful asymptotics. These review and extend a variety of 
results already in the literature. 

For definiteness, we consider pulsars with an acute angle x be- 
tween S2 and B (‘acute’ pulsars). These objects have < 0 at the 
polar cap, and require electron emission to supply jjgj- 

Let us consider an electron beam starting with zero velocity at 
x = 0 and let the current density j m imposed by the magnetosphere 
(equation (2)) be a fraction f of the GJ charge density 

jm = ? Jgj • (5) 

In our case the GJ charge density is negative. In stationary flow the 
current density is constant in both space and time and is equal to the 
imposed current density j = j m , thus dE/dt = 0 in equation ( 1 ). The 
stationary electric field E s is then given by Gauss’s law which in 
the frame corotating with the NS takes the form (e.g. Fawley et al. 
1977, and references therein) 

d E s 

— = 4n(r] - jj G j) ■ (6) 

ax 

The magnitude of the electric field increases with distance if the 
magnitude of the charge density rj is larger than the GJ charge 
density and decreases otherwise. 


2 TV potential drops are required if curvature radiation from charges ac- 
celerated along a locally dipole magnetic field is the emission mechanism 
for the gamma-rays that convert to pairs. The required potential drops are 
smaller if the emission mechanism is inverse Compton scattering (either res- 
onant or non-resonant) of softer photons from the surface (e.g. Hibschman 
& Arons 2001, and references therein). 
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The charge density at any given point of the flow is 


I = j/v = 


jm V p 2 + 1 
c p 


(7) 


where v is the flow velocity and p = y v/c is the four- velocity = 
momentum in units of me. p of a charge-separated stationary beam 
is given by the solution of the equation (BIO) 


( 



= 2 P + 1 (l+%p- vV + l) , 


( 8 ) 


where the distance s is measured in units of the Debye length k Di gj 
of a cold electron plasma with GJ number density 


7-d.gj 


c ^ 47D; GJ e ^ 


-1/2 


~2 B,T 1/2 P 1/2 cm. 


(9) 


B 12 is the pulsar magnetic field in 1() 12 G, and P is the pulsar period 
in seconds. 

The numerical solutions of (8) for different imposed current den- 
sities are shown in Figs 2 and 3 for different values of these 
results confirm the work of Shibata (1997). For 0 < j m /jo] < 1 the 
steady flow oscillates spatially, with particle momenta oscillating 
in the interval [0, p max ], with 

2 ? 

Pmax = Tn- (10) 

1 - 

dp/ds = 0 at p = p max and so is the RF1S of equation (8). The 
value of p max and the spatial period of oscillations jo increase with 
increasing £ (see Figs 2 and 3). For § > 1 acceleration is monotonic 
with p increasing to infinity, with asymptotic behaviour 

p = Vis + ^-^-s 2 . (11) 

See Fig. 3 for£ = 1, 1.1, 1.5. 

The reason for such behaviour is as follows. The flow starts 
with zero initial velocity at x = 0 where the electric field is zero. 



Figure 3. Phase space trajectories for stationary space charge limited flow 
for current densities jm/jai = 0.9, 0.95, 0.99, 1, 1.1, 1.5. Note that in contrast 
to Fig. 2 the vertical axis on this plot is logarithmic. 

When the particle velocity is small the imposed current density is 
produced by high particle density moving slowly. In such places 
the absolute value of the beam charge density (cf. equation 7) is 
larger than that of the GJ charge density, |t/| > |?/gjI (see Fig. 4 
where we plot the ratio for flows with § = 0.5, 1, 1.5). The 
charge density is negative and according to equation (6) the electric 
field in this region is decreasing towards more negative values, thus 
accelerating electrons. If the imposed current density exceeds the 
GJ current density, f > 1. the absolute value of the beam charge 
density - whose maximum value is [/ m /c| = £ I 7 gj I - never becomes 




Figure 2. Phase space trajectories (four-velocity p versus distance normal- 
ized to /.[) G j) for stationary space charge limited flow for current densities 
jm/ /GJ = 6.1, 0.25, 0.5, 0.75, 0.9, 0.95. 


Figure 4. Charge density of stationary space charge limited flow normalized 
to the GJ charge density r]Q j as a function of distance x normalized to /,□ 
for current densities jm/jai = 0.5, 1.0, 1.5. 
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smaller than [ rjc j | , hence, dE s /dx < 0 and acceleration continues up 
to infinity- the electric held and potential are monotonic. For f < 1, 
on the other hand, particles’ velocities increase to values such that 
the imposed current density can be sustained by particle number 
density smaller than the GJ number density. Then \ r]\ < [jjgjI ( see 
the line for = 0 . 5 /gj in Fig. 4 ), dE s /dx > 0 and the accelerating 
electric field weakens, changes sign and decelerates particles - in 
cold flow, the particles decelerate to zero velocity, and the cycle 
repeats. 

The model of space charge limited flow outlined here provides a 
physical framework for the expected particle energetics. It is based 
on the approximation of one cold fluid and an assumption of com- 
plete stationarity of the flow. It can be extended to a two-fluid 
model to account for the presence of positrons and pair creation (as 
in Arons 1983). However, kinetic effects such as particle trapping 
cannot be included in a cold fluid approximation, although certain 
aspects can be modelled if momentum dispersion (‘pressure’) is 
included, with an assumed equation of state. A kinetic model incor- 
porates momentum dispersion in the collisionless medium without 
having to make arbitrary assumptions about the equation of state. 
As we show in the following sections, particle trapping and pair 
creation profoundly affect the plasma dynamics, with momentum 
dispersion being essential to the dynamics behind the simultaneous 
adjustment of the charge density to the condition of low voltage drop 
along B, modelled as E ■ B = 0 in the force-free global model, and 
the adjustment of the field aligned current j to the magnetospheri- 
cally imposed j m . 

The study of plasma kinetics in a general regime - without relying 
on stationarity or perturbation theory - is possible only by means 
of numerical simulations. In following sections we describe our 
study of plasma, both fully non-neutral and quasi-neutral when pair 
cascades form, with the help of a self-consistent hybrid numerical 
model incorporating both charged particles and photons. 


4 NUMERICAL SETUP 

We use the same one-dimensional hybrid Particle-In-Cell/Monte 
Carlo hybrid code described in Timokhin (2010) modified for the 
space charge limited flow regime. Below we briefly describe the 
main equations, notations and numerical algorithms; a detailed de- 
scription can be found in Timokhin (2010, sections 2 and 3). 

We solve the evolutionary equation for the electric field E 

9 E(x, t ) 

a = —471 OX*, D ~ jm) , (12) 

at 

where j(x, t) is the actual current density and j m is the current 
density imposed on the cascade zone by the magnetosphere. This 
equation is Ampere’s law, equation (1). We are solving an initial 
value problem; thus an initial distribution of the electric field E{x, 
t = 0) must be supplied. At the start of the simulation we construct 
the initial distribution of the electric field by solving the Gauss 
equation for the electric potential 4> assuming some initial charge 
density distribution r/ start 


d 2 (p 

dx 2 

E = 


-47t(>J star t - t] Gj) 

dtp 

dx 


(13) 

(14) 


We proceed with the time integration of equation (12) using a charge 
conserving scheme (e.g. Birdsall & Langdon 1985; Villasenor & 
Buneman 1992), so the Gauss equation is satisfied at each successive 


time-step up to machine precision. The GJ charge density enters in 
equation (13) for the initial configuration of the electric field; this 
information is then ‘carried on’ in time by equation (12). 

To model the space charge limited flow at every time-step we 
inject electrons and protons just outside the numerical domain used 
for the electric field calculation and let the system pull the necessary 
amount of particles into that domain. We do not set E(x = 0, t) to 
zero as a boundary condition but rather allow the plasma in the 
system to enforce this condition as part of the simulated physics. 
A detailed description of our algorithm for reproducing the space 
charge limited flow condition at the NS surface is given in Appendix 
C. When pair creation cascades occur, we take into account only 
curvature radiation as the gamma-ray emission mechanism; pairs 
are created by single photon absorption in the strong magnetic field 
(e.g. Erber 1966). 

We performed many numerical experiments starting from differ- 
ent initial conditions: (i) computational domain filled by plasma 
with charge density equal, less and higher than the GJ charge den- 
sity as well as starting with vacuum; (ii) different initial potential 
drop over the domain; (iii) different length of the computation do- 
main. In all cases without exception after initial relaxation on the 
time-scale of the order of the fly-by time of the domain the system 
settled down to a configuration which depends only on the imposed 
current density j m . 

5 LOW ENERGY CHARGE SEPARATED, 

SUB-GJ FLOW: 0 < j m /j GJ < 1 

As it turns out, there is no pair formation for 0 < j m /joj < 1 
and the only characteristic spatial scale of the flow is the Debye 
length kn, gj- In this section we will discuss the properties of such 
flow using simulations in domains with the length L up to several 
hundreds of k D , gj, which is much less than the width of the polar 
cap. Using such a small domain we well resolve the characteristic 
spatial scale; increasing the domain length does not change the 
results. 

According to the stationary solution from Section 3 there is no 
relativistic particle acceleration when 0 < j m /jo] < L The flow is 
spatially oscillatory and the maximum momentum of particles p max 
is by far too small for emission of pair producing photons. It is 
unlikely that such oscillatory cold flow can exist - near the stagna- 
tion points it may be a subject to ‘instability’ with characteristics 
of the wave breaking to which non-linear waves in cold fluids with 
velocity stagnation are subject, which could destroy the spatial os- 
cillation and create an effective resistor, across which a substantial 
fraction of the perpendicular (to B) voltage drop might appear. The 
available voltage drop across B is huge and in principle it might 
happen that the system ends up in a state with highly oscillatory 
electric field and bursts of pair formation as predicted in the two- 
fluid plasma model of Levinson et al. (2005). On the other hand, 
finite amplitude electrostatic waves containing similar stagnation 
points show wave breaking, with alteration of cold to hot flow with 
momentum dispersion no more than comparable to the flow veloci- 
ties found in the originally constructed cold oscillation (Akhiezer & 
Polovin 1956; Tajima & Dawson 1979). That would create a warm 
but non-relativistic or mildly relativistic charge-separated outflow. 

We find that in the sub-GJ regime, the non-neutral flow is indeed 
low-energetic, with particle energies orders of magnitude below 
the energy/particle required for pair production. However, the final 
state differs drastically from the oscillatory flow from the stationary 
solution shown in Figs 2 and 3. Even if at the beginning of the 
simulations particles follow trajectories of the stationary solution, 
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Figure 5. Development of ‘trapped particle’ flow when the space charge limited flow which starts into vacuum. Phase space portrait of particles are shown for 
16 moments of time indicated in small boxes on the top of each plot. The current density j m = 0.5_/gj- The distance x is measured in units of the Debye length. 
Red dashed lines show the analytical solution for stationary flow. Particle momenta p_ are normalized to m e c. The total length of the computational domain 
L = 50/d, gj, only part of it ( x < 50A.d, gj) is shown here. Time is measured in fly-by time (L/c) of the whole domain. Snapshots in the same row have the 
same time interval between them, but these time intervals are different for different rows; they increase towards the bottom row. 


the standing non-linear wave structure breaks quickly - particles at 
the velocity zeros of the cold flow go both up and down. 

A good example of this inherent instability is shown in Fig. 5. 
We start from a vacuum configuration - when there are no particles 
in the domain - and let the system evolve. In Fig. 5 we show 
snapshots of the phase space portraits of the flow with j m = 0.5/gj; 
the distance from the NS, normalized to ko, gj, is along x-axis, and 
particle momenta, normalized to m e c, are along y-axis. The whole 
domain has the length L = 50ko, gj an d we show only a part of 
it here. Time t is measured in fly-by time of the whole domain 
L/c. With the dashed red line we show phase space trajectories of 
particles from stationary solution. Particles coming from the surface 
at first follow the trajectories of the stationary solution. However, 
after coming to the first stagnation points some of the particles are 
turned back and the flow starts to randomize. After several tens of 
plasma periods k D , gj /c the flow reaches its final configuration. 

Examples of final configurations for space charge limited flow 
with different current densities are shown in Fig. 6, where we plot 
phase space portrait of particles in the whole computation domain. 
The flow has two components: a warm beam of particles with high- 
est momentum which produce the required current density and a 
cloud of charged particles circulating in the domain - these com- 
pose an electrically trapped, ‘thermal’ component. In the cloud 
component there is roughly equal number of particles moving in 
opposite directions; these particles do not contribute to the current 


but contribute r ] gj — j m /c to the charge density keeping t/ tota i equal 
to the GJ charge density. The distinction between these components 
is not absolute as some particles from the beam go into the cloud and 
vice versa, although the fraction of mixing particles is small. Some 
of the particles in the cloud have very low momenta and so they can 
adjust to any given charge density. Hence, in sub-GJ space charge 
limited flow, 0 < j m / y/jj < 1, the electric field is not sensitive to vari- 
ation in the GJ charge density 3 - in contrast to the large importance 
of variation in the GJ charge density for relativistic acceleration of 
space charge limited cold beams in the polar cap cascade models of 
Arons & Scharlemann (1979) and Muslimov & Tsygan (1992). 

Plasma flow in the sub-GJ regime can be described as a beam of 
mildly relativistic particles propagating through a cloud of trapped 
particles with near-thermal distribution. In Fig. 7 we plot particle 
distribution functions. The beam component is visible as a bump 
on the distribution function at the high momentum side. The cloud 
component has a near-thermal (Maxwell-Juttner) momentum dis- 
tribution (at least in its low-energy part) dn/dp oc const - such 
quasi-thermalization is a common consequence of the phase mix- 
ing between the particles and fields built into the wave breaking 
process. 

3 We performed simulations for jm/jai < 1 with variable GJ charge density, 
and, as expected, saw the electric field to be just as screened as in the uniform 
GJ density case. 
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Figure 6. Phase space portraits of well-developed space charge limited flows for six different current densities: jm/jai = 0.1, 0.25, 0.5, 0.75, 0.9, 0.95. Current 
density j m is indicated in the left upper part on each plot. Distance is measured in units of a\). and particle momenta p are normalized to m e c. Note that the 
lengths of the computation domain differ between these plots. 



Figure 7. Particle distribution functions p ( dn/dp ) for well-developed space charge limited flows from Fig. 6. Distribution functions of particles with positive 
momenta (moving towards the magnetosphere) are shown by solid lines, and distribution functions of particles with negative momenta (moving towards the 
NS) by dashed lines. 


When particles leave the NS they are non-relativistic and their 
charge density p = j m /v = ^joi/ v is larger than the GJ charge 
density and so they form a charge sheet near the surface generating 
accelerating electric field, just as in the idealized stationary case (see 
Fig. 4). When particles reach the velocity such that \j m \/v < [??gjI 
the electric field derivative dE/ds changes sign and after some dis- 
tance the electric field can start decelerating particles. At the point 
where E = 0 particles reach their maximum momentum. Above the 
gap particles from the cloud component add additional charge and 
so the change density there is equal to jjgj and the electric field is 
screened. The length of this gap is of the order of .v max = so/2 (half 
the spatial period of cold flow oscillations) and so the maximum 
momentum particles gain in this gap is comparable to p max (both 


■''max an d Pmax depend on j m ). In Fig. 9 we plot momenta of particles 
in the beam component as functions of § = j m /joi superimposed 
on the theoretical dependence of p max given by equation (10); the 
agreement is pretty good. The current density in the final configu- 
ration is close to j m throughout the domain, with small fluctuations 
around this value Sj <g. j m . 

In Fig. 8 we plot the electric field in the calculation domain at the 
same moments of time as the phase space portraits shown in Fig. 6. 
The electric field at any given point fluctuates, but the relative fluc- 
tuation at the beginning and at the end of the domain (accelerating 
and decelerating regions; see below) are much smaller than those in 
the centre of the domain (the region of the charge cloud). The elec- 
tric field in the gap near the NS launches the beam component. The 
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Figure 8. Snapshots of the electric field E for well-developed space charge limited flows taken at the same moment as the phase space portraits from Fig. 6. 
E is normalized to 7T?7gj^d, gj- 



Figure 9. Momenta of the current currying particle beam. The widths of 
the peaks at Fig. 7 as the function of f = j m /jas arc shown by the blue bars. 
The solid red line is the maximum achieved momentum for the stationary 
SCLF solution. 

electric field at the other end of the computation domain supports 
the cloud components by reversing the momenta of most of the par- 
ticles in the cloud moving away from the NS, sending them back. 
This electric field is not strong enough to reflect most of the beam 
particles. The mixing between the beam and cloud components is 
the strongest here. This electric field appears self-consistently when 
the flow reaches its steady configuration, 4 because it is needed to 
sustain the cloud component necessary to match both the charge and 
the current density in the domain. In our small-scale ID simulations 
we impose the current density on the domain, which finally gives 
rise to the electric field at the right end of the domain. In reality, it 
is the magnetosphere which sets the current density by twisting the 
magnetic field lines and generating electric field reversing some of 
the particles. This second region with an unscreened electric field 
at the magnetosphere end of the domain in our model may be a 

4 Th e formation of the cloud component is not linked to the appearance of 
the electric field at the end of the domain (see e.g. Fig. 5). 


‘compressed’ version of some parts of the outer magnetosphere. 
For example, on field lines that pass through the null surface where 
Si ■ B = 0, the plasma cloud and beam, composed of only one sign 
of charge, cannot freely enter the outer magnetosphere (Scharle- 
mann, Arons &Fawley 1978) in the absence of other sources of elec- 
tric field in other parts of the magnetosphere (Goldreich & Julian’s 
‘hanging charge clouds’), offering the possibility of opening a vac- 
uum gap within the otherwise force-free structure. On the other 
hand, on polar field lines that never cross the null surface - most of 
them, in the aligned rotator - the charge-separated beam and cloud 
can extend outwards ‘forever’, in principle. Multi-dimensional par- 
ticle simulations, analogous to those of Spitkovsky & Arons (2002), 
are required to see if indeed this speculation is true (as well as de- 
termine how this essentially ID model might fit together with the 
other, more ‘lively’ aspects of the polar flow outlined in Section 6). 

In Fig. 10 we show power spectra of the fluctuating electric field 
in the central parts of the calculation domain, outside of the accel- 
eration zones, p = \Eh\ 2 , where is the spatial Fourier amplitude 
of the electric field and the wave vector k is normalized to the Aqqj. 
These spectra can be fitted with the power law p oc kr a with a 
between 2 and 3 for all current densities. 

Thus, along magnetic field lines where the current density is sub- 
GJ there is no pair formation, so long as strong electric fields in 
neighbouring, more active flow zones do not leak into the cloud. In 
an aligned rotator, for example, no pairs are forming above most of 
the polar cap area. Plasma flowing along such magnetic field lines is 
mildly relativistic, consists of particles of only one sign (electrons) 
and its density is low, equal to the GJ number density. However, 
observational evidence for ongoing pair formation in pulsars is very 
strong and, hence, there must be regions in the pulsar magnetosphere 
where electron-positron plasma is created. 

6 PLASMA FLOW WITH PAIR FORMATION - 
SUPER-GJ AND RETURN CURRENT REGIONS 

Let us now consider what happens along magnetic field lines where 
the current density has either (i) the opposite sign to the GJ current 
density, j/ ja < 0, anti-GJ flow, or (ii) the same sign and absolute 
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Figure 10. Power spectra of the electric field /* = |£r | 2 for well developed space charge limited flows from Fig. 6. k is normalized to the corresponding 
IAd, gj- 


value larger than the GJ current density, j/jcs > 1, super-GJ flow. 
Case (i) includes the regions of return current, including the current 
sheet, while case (ii) is relevant for most of the magnetic field lines 
in a nearly orthogonal rotator. As we show in this section, in both of 
these cases a strong accelerating electric field is generated. When 
the resulting potential drop is sufficient to accelerate particles up to 
the energies such that they can emit photons that convert to pairs 
within the accelerator, the plasma flow will be highly non-stationary 
with intermittent pair creation. 

In a real pulsar the magnetic field strength falls with the distance 
(oc r~ 3 for a dipole field) and pair creation is possible only at 
sufficiently low altitude (if gamma-ray interaction with thermal low 
energy photons is neglected). In order to imitate the effect of pair 
creation attenuation with the distance we set the magnetic field 
strength to zero starting at distance x B from the NS - the magnetic 
field is given by 


B(x) = 


So, 

0 , 


if x < x B 
if x > x B . 


(15) 


If the charge density were formed from steadily flowing rela- 
tivistic beams, it would vary with the height as rj(x) oc B(x). The GJ 
charge density changes in a slightly different way, r](x) oc B(x)f(x) 
due to inertial frame dragging (Muslimov & Tsygan 1992; Beskin 
1990) or/and field line curvature (Scharlemann et al. 1978). If one 
neglects the latter effects, scaling oc B(x) can be incorporated into 
the spatial (for equation 13) or temporal (for equation 12) coordi- 
nates. The electrodynamics of the cascade zone can be modelled in 
a ID problem with constant GJ charge density; the only effect of the 
GJ charge density variation will be in changing the spatial (and tem- 
poral) scales. The same scaled ID model can also be used to study 
the effects on the electrodynamics of the cascade zone of deviation 
of the GJ charge density scaling from being oc B{x) by consider- 
ing a problem where gj depends on the distance; the variation of 
the GJ charge density will be given by f(x) as the dependence on 
B(x) is already incorporated in the model. First, in Sections 6. 1 and 
6.2, we consider the case when the GJ charge density is constant, 
i.e. this case corresponds to a model where variation of rjcj/B is 


neglected. Then in Section 6.3 we address the influence of the GJ 
charge density variation on the physics of the cascade zone. 

While simulations of space charge limited flow with sub-GJ cur- 
rent densities described in Section 5 can be considered as directly 
related to more complete pulsar models - the distance over which 
wave breaking and trapping control the flow is much smaller than 
the width of the polar cap, and the ID approximation is well mo- 
tivated - the pair creation models presented in this section can be 
considered only as illustrative of the physics but not fully applica- 
ble to pulsars. The spatial scales over which pair creating photons 
are absorbed, even when they are emitted at very low altitude (e.g. 
the attenuation length of the magnetic field), are much larger than 
the width of the polar cap, making transverse structure essential for 
modelling the pulsar environment - such effects are necessary for 
a full evaluation of the pair yield, since much of the pair creation 
occurs in regions beyond the acceleration zone. Nevertheless, from 
our ID models we obtain insight into the basic cascades physics 
in a regime never modelled previously. Multi-dimensional models 
will be considered elsewhere. 

Within the context of the ID model, the domain length, L, and 
the value of x B (the height above which the magnetic field is too 
weak to support pair creation) are the parameters having the largest 
departure from what would appear in multi-dimensional context. 
These lengths in our simulations are much smaller than in real 
pulsars, but that choice allowed us to construct models utilizing 
reasonable amounts of CPU time and to explore a broad parameter 
space. L and x B were chosen in such a way that the minimum size 
of the accelerating region necessary to start pair creation is at least 
two to three times less than .v B . 

We performed numerous simulations for different initial con- 
ditions and physical parameters in order to study the qualitative 
behaviour of cascades. We performed simulations with different 
values of the numerical parameters (spatial resolution, time-steps, 
particle injection rate, number of particles per cell) in order to 
check the numerical convergence. In all physical cases presented 
in this section plasma flow is quasi-periodic, and this behaviour 
does not depend on initial conditions - after a short relaxation time, 
comparable to the fly-by time of a relativistic particle through the 
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computational domain, the system settles down to a limit cycle 
sort of behaviour. We describe here a particular set of simulations 
which is representative of all other models. Simulations described 
in Sections 6.1 and 6.2 have the following physical parameters: 
the length of the domain L = 2.4 x 10 4 cm, the potential drop 
in vacuum across the domain AV = 10 14 V, 5 the radius of curva- 
ture of the magnetic field lines with respect to a photon orbit p c = 
10 6 * * * * * cm ~ the stellar radius (small compared to the pure star cen- 
tred dipole value ~/R*R lc ~ 1 0 7 S P 1/2 cm but easily attained in 
offset dipole geometry (Arons 1998; Harding & Muslimov 2011; 
Arons, in preparation); magnetic field strength So = 10 12 G. The 
distance x B marking the transition to the outer magnetosphere with 
a small magnetic field is set to x B = 0.7 L. The particular simulations 
described in these two sections differ only in the imposed current 
density j m . 

We illustrate the flow dynamics with a series of snapshots for 
different physical quantities shown in Figs 11-16 and 20-22. In 
these figures, each column gives detailed information about physical 
conditions in the computation domain at a given moment of time: the 
number densities of electrons and positrons n± (shown as charge 
densities of electrons and positrons p±, n± = |^±|), total charge 
density p, current density j, the accelerating electric field E. phase 
portraits of electrons, positrons and pair producing photons, and, 
in Figs 14-16, protons. Particles with positive values of the four- 
velocity p are those which move away from the NS (towards higher 
altitude), and particles with negative/; move towards the NS. These 
plots are similar to ones in Timokhin (20 10), with the only difference 
that now we use semi-logarithmic scale for particle momenta (linear 
for —5 < p < 5 and logarithmic everywhere else) on phase space 
portraits that show dynamics of high and low energy particles on 
the same plot. 

The number density, charge density and the current density are 
normalized to the corresponding GJ values: p± and p are normal- 
ized to I Pqj \ , j to |)/gj|c, where Pqj is the GJ charge density at the 
NS surface (the distinction between p gj and Pq S will be important 
for models described in Sections 6.3 and 6.4). The electric field 
is normalized to Eg = \Pqj\uL. The distance x on these plots is 
normalized to L, much larger than the Debye length A. Di gj- The 
time t is normalized to the relativistic fly-by time of the computa- 
tional domain L/c = 0.8 p.s for the chosen parameters. The time 
is counted from the start of a particular simulation, so only time 
intervals between the snapshots have a physical meaning. 

In all cases, plasma flow and pair formation have limit cycle 
behaviour. In each case we illustrate this behaviour by three series 
of snapshots taken within one typical cycle. These three series show 
three phases of plasma flow: cascade ignition (Figs 11, 14 and 20), 
development of the cascade (Figs 12, 15 and 21) and filling the 
domain with dense pair plasma (Figs 13, 16 and 22). In each figure 
time intervals between snapshots are equal, but these intervals are 
different in different figures. 


5 Such vacuum potential drops over the domains of this size are realistic 

in young, high magnetospheric voltage pulsars, O nl > 10 15 V. This choice 

allows studying pair formation over distances smaller compared to the polar 

flux time diameter. In more common pulsars with < 10 13 V, pair for- 

mation happens over (much) larger distances. However, as we mentioned 

before, our simulations represent a toy model addressing the general be- 

haviour of the polar cap acceleration zone and our choice of parameters is 

motivated by convenience of simulations, rather than by attempt to model a 

real pulsar. 


6.1 Flow with y m // GJ < 0 

In this section we describe cascade development for two cases when 
the imposed current density is anti-GJ (i.e. has the opposite sign to 
the GJ current density), for j m = — 0 . 5 / 0 ! and jm = — 1 - 5 / gj . To sup- 
port such current electrons on average must move towards smaller 
x, down towards the NS, positrons towards larger x, up into the mag- 
netosphere. Electrons could be freely emitted from the NS surface, 
but because the imposed current causes the growing electric field 
to point away from the surface, electrons at the surface are acceler- 
ated downwards and none is extracted from the star. 6 This situation 
resembles the case of the Ruderman & Sutherland (1975) cascades 
studied in Timokhin (2010) in the sense that all the particles are 
produced during the burst of pair formation. Hence, the physics of 
pair cascades with anti-GJ imposed current density discussed in this 
section are also applicable to Ruderman & Sutherland cascades as 
well. 

Particles leave the domain and some time after the burst of pair 
formation the domain becomes depleted of electrons. This process 
starts at the right end of the domain - at the ‘magnetosphere’ end 
- as electrons are moving down towards the NS. In the region 
depleted of electrons, the positrons support the current density j < 
j m with charge density p > paj. This gives rise to the electric field 
which accelerates positrons towards the magnetosphere - see the 
phase portraits of positrons in Figs 11 and 14. The electric field 
in this region at any given point is growing with time and the size 
of the region with unscreened electric field is getting bigger as the 
remaining electrons are moving towards the NS. The electric field 
grows linearly with the distance because positrons are relativistic 
and so their charge density remains constant (see plots for E in Figs 
1 1 and 14). 

As the region with the unscreened electric field grows, positrons 
are being accelerated up to higher and higher energies and begin 
emitting pair producing gamma-rays. Positrons remaining from the 
previous burst of pair formation are the particles which ignite the 
next burst. As positrons are moving up (and so do the first pair pro- 
ducing photons), the first pairs are produced at the largest distance 
from the NS where pair formation is still possible, in our case near 
x = x B = 0.7L. The injected pairs are picked up by a very strong 
electric field and are accelerated to high energies in less time than 
the first generation positrons. They emit pair-producing capable 
gamma-rays, but now both secondary electrons and positrons are 
emitting pairs. As electrons and positrons move in opposite direc- 
tions so do the gamma-rays - see snapshots for gamma-rays phase 
space at f > 8.100 in Fig. 11 and t > 6.390 in Fig. 14. 

Secondary electrons and positrons are moving in opposite di- 
rections and the plasma gets polarized - see plots for p± and p 
at t = 8.260 in Fig. 11 and at t = 6.430 in Fig. 14. When their 
number density becomes comparable to the GJ number density, 
these particles start screening the electric field (see plots for E at 
t = 8.260 and 8.270 in Figs 1 1 and 12 and at t = 6.430 and 6.440 
in Figs 14 and 15). The screening starts in the region where the first 
pairs were injected because pairs have been injected here for the 
longest time and their number density is larger. 

The rate at which particles left from the previous burst of pair 
formation are leaving the domain depends on the imposed current 
density. The average bulk motion of electrons for j m = — 0 . 5 /gj is 

6 Due to numerical noise the electric field fluctuates and sometimes elec- 
trons ‘from the NS surface’ enter the domain; however, the number of such 
electrons is well below the fluctuation of electron number density due to 
numerical noise. 
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Figure 11 . Ignition of pair formation in anti-GJ flow with j m = — 0.5/gj- Several physical quantities are shown as functions of the distance x from the NS; 
x is normalized to the domain size L. Plots in each column (for the same time f) are aligned - they share the same values of x. The following quantities are 
plotted: first row: p\ — charge density of electrons (negative values, blue line) and positrons (positive values, red line); rj± is normalized to the absolute value 
of the Goldreich-Julian charge density | r?GJ I - Second row: the total charge density p normalized to the absolute value of the Goldreich .lulian charge density 
|)7 gj|. Third row: current density j normalized to the absolute value of the Goldreich-Julian current density [/gj I = |t?GJc|. Fourth row: accelerating electric 
field E normalized to the ‘vacuum’ electric field Eq = | (Jgj 1 7t L. Fifth row: phase space portrait of positrons (horizontal axis - positron position x, vertical axis 
- positron momentum p \ normalized to m e c). The vertical axis is logarithmic except for the region around zero momentum (—5 < p \- < 5), where the scale is 
linear. Sixth row: phase space portrait of electrons (horizontal axis - electron position x, vertical axis - electron momentum p normalized to m e c). Seventh 
row: phase space portrait of pair-producing photons (horizontal axis - photon position x , vertical axis - photon momentum p y normalized to m e c). 


less than 0.5c, while for j m = — 1.5/gj it is close to c, 7 and the size 
of the electron-depleted region grows slower in the first case than in 
the second one. Second generation electrons, which mark the upper 
boundary of the gap, move with ultrarelativistic velocity towards 
the NS. Because of this in the case of j m = — 0.5/gj the gap with the 
accelerating electric field quickly disappears - at t= 8.430 in Fig. 12 


7 Note that the time intervals between the snapshots in Fig. 1 1 is two times 
larger than that in Fig. 14. 


secondary particles have already caught up with the particles from 
the previous burst of pair formation and the electric field is screened 
everywhere where pair formation is possible. In the case of j m = 
— 1.5/gj the gap moves towards the NS approximately retaining its 
size (see Figs 15 and 16). This behaviour is generic. Namely, for 
imposed current densities with less-than-GJ absolute values, [/ m | < 
[/"gj I , the average bulk motion of particles from the previous burst 
of pair formation is non-relativistic, which leads to quick closure 
of the accelerating gap. If the absolute value of the imposed current 
density is larger-than-GJ, [/ m | > [/gj I , the average bulk motion of 
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Figure 12. Screening of the electric field in anti-GJ flow with j m = — 0.5/gj- The same quantities are plotted as in Fig. 11. 
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particles from the previous burst of pair formation is relativistic and 
the gap propagates in the direction of this bulk motion. The same 
behaviour was observed also in the case of Ruderman-Sutherland 
cascades studied in Timokhin (2010). 

In the case of j m = — 0.5/gj the gap with the accelerating electric 
field closes and particles in the region where pair formation is possi- 
ble are not accelerated anymore. The pair formation continues only 
due to particles being accelerated before the gap disappeared. There 
are comparable numbers of particles accelerated in both directions, 
but electrons moving towards the NS propagate in the region where 
pair formation is possible and so after screening of the electric field 
most of the pairs are created with the initial momenta directed to- 
wards the NS (see phase portraits for pair producing gamma-rays 
for t > 8.180 in Figs 11-13 and particle distribution functions in 
Fig. 17). In the case of j m = — 1.5/gj the gap propagates down 
to the NS surface. While the gap is moving, positrons left from 
the previous bursts of pair formation, which are still present at 
the NS side of the domain, enter the gap, are picked up by the 


electric field, and emit pair producing gamma-rays. In this case 
the pair production is sustained not only by electrons accelerated 
in the initial screening event but also by positrons continuously 
accelerated in the moving gap towards the magnetosphere. Phase 
portraits of gamma-rays in Figs 14-16 clearly show the presence of 
gamma-rays moving towards the magnetosphere during all stages 
of cascade development (see also the plot for particle distribution 
functions in Fig. 18). Gamma-rays with positive momenta occupy a 
larger fraction of the space over time, as the gap propagates towards 
the NS. 

Regions far from the NS surface experience charge starvation 
first. The pair formation is possible only in the regions with a strong 
magnetic field and until the starved region extends into the strong 
field domain the electric field remains unscreened. After the start of 
pair formation, the dense plasma propagates also in the direction of 
the magnetosphere screening the accelerating electric field there. In 
our simulations the region with x > 0.7 L represents the rest of the 
magnetosphere, and in snapshots with t = 8.430-8.990, in Figs 12 
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Figure 13. Filling of computation domain with dense pair plasma in anti-GJ flow with j m = — 0.5/gj. The same quantities are plotted as in Fig. 11. 


and 13, and t = 6.560-6.930, in Figs 15 and 16, one can clearly 
see how the dense pair plasma fills the domain with x > 0.7 L. 
Filling of that domain with plasma is forced by a small electric field 
induced in the dense pair plasma by the imposed current. The rate 
at which these regions are filled with the dense plasma depends on 
the imposed current - in the case of j m = — 0.5/gj the average bulk 
motion of the pair plasma is about 0.5c, while for j m = — 1 ,5jo j it is 
relativistic. 

The accelerating electric field is very strong in the domain where 
pairs cannot be produced until electrons generated in the discharge 
arrive, and positrons (either primary or secondary ones) entering this 
domain before electron arrival are accelerated up to high energies. 
In our model pairs are created only by single photon absorption in a 
strong magnetic field, which automatically restricts pair formation 
to the polar cap region. However, electron-positron creation via 
y-y absorption can be possible at much higher distances from the 
polar cap, similar to the outer gap scenario (Cheng et al. 1986). 


In our model between the bursts of polar cap pair formation in 
the magnetospheric region above the low-altitude pair producing 
zone, atx> Xb, particles can be accelerated up to radiation-reaction 
limited energies and emit very high energy gamma-rays. In a more 
realistic model, when y-y absorption is taken into account, this 
should give rise to cascades in the outer magnetosphere resembling 
to some extent the outer gap scenario, except that the accelerating 
zone would not be limited to the region around the null surface, 
where the GJ charge density changes sign. Due to the intermittency 
of plasma flow in the return current region both polar cap and the 
outer gaps like cascades might exist along the same magnetic field 
lines. However, without higher dimensional simulations it is not 
clear how the interaction between two such cascade zones would 
play out. 

Screening of the accelerating electric field during each burst of 
pair formation gives rise to a superluminal electrostatic wave. In 
Fig. 19 we show screening of the electric field in the gap for the 
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Figure 14. Ignition of pair formation in anti-GJ flow with j m = — 1 .5/gj . The same quantities are plotted as in Fig. 1 1 with the addition of the phase space 
portraits for protons in eighth row: (horizontal axis - proton position x , vertical axis — proton momentum p p normalized to ;n p c). 


flow with j m = — 0.5/gj, for flow with j m = — I .5/gj the situation is 
very similar. Plots in Fig. 19 show in more detail than in Fig. 12 how 
the accelerating electric field is being screened by newly created pair 
plasma about t = 8.350. Screening of the electric field starts in the 
middle of the blob of newly created plasma and spreads to its edges. 
This spreading occurs in the form of an electrostatic wave the phase 
velocity of which is larger than c. This process is almost identical 
to what was observed for cascades described in Timokhin (2010) 
where more detailed discussion of the screening process can be 
found in section 4.2 [cf. our Fig. 19 with fig. 5 in Timokhin (2010)]. 

As in the simulations described in Timokhin (2010), second gen- 
eration particle momentum distributions are very broad, extending 


towards non-relativistic energies. At least one of the reasons for 
momentum broadening of the particle spectra is trapping of parti- 
cles in the strong electrostatic wave excited at the beginning of each 
burst of pair formation. The dynamics of these discharges are very 
similar to what is seen in simulations of the Ruderman & Suther- 
land cascade; more details on momentum spreading can be found in 
section 4.3 of Timokhin (2010). The appearance of the low energy 
component in the pair plasma is visible in the electron and positron 
phase portraits - initially at distances where pairs are injected, there 
are no particles in the phase space with low momenta, but later 
particles fill in the low momentum regions (see plots at t > 8.350 in 
Figs 12 and 13 and for t > 6.500 in Figs 15 and 16). In Figs 17 and 
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Figure IS. Screening of the electric field in anti-GJ flow with ;' m = — 1.5/gj- The same quantities are plotted as in Fig. 14. 


18 we plot the particle momentum distribution p (dn /dp) for three 
different moments of time. In the upper panel we plot the momen- 
tum distribution of particles moving towards the magnetosphere, p 
is positive; in the lower panel - the momentum distribution of par- 
ticles moving towards the NS, p is negative. These distributions are 
averages for regions where the electric field at the magnetosphere 
end is screened. On these plots one can see that the low energy 
plasma component is present at all stages of cascade development. 

There is a physical effect specific to cascades along field lines 
carrying return current - positive ions (in our model protons) can 
be pulled from the NS. Both electrons and protons can be extracted 
from the NS surface, and in our simulations we allow injection 
of both species. The imposed current density requires positively 


charged particles to move towards the magnetosphere, i.e. protons 
could be pulled from the NS. However, for less-than-GJ current den- 
sities ([/ m | < L/gjI) there are always electrons and positrons near the 
NS surface left from the previous burst of pair formation. The gap’s 
accelerating electric field disappears before it reaches the NS, while 
pair plasma appears after every burst of pair formation keeping the 
electric field near NS screened. In our simulations of cascades with 
less-than-GJ current densities we did not see injection of protons. 
In the case of larger-than-GJ current density ([/ m | > L/'gjI) protons 
are indeed pulled from the NS, at least at some stage of cascade de- 
velopment, and this cannot be attributed only to numerical effects. 
In this case we see some protons pulled from the NS at any time, 
even if pair plasma is present near the NS. The number density of 
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Figure 16. Filling of computation domain with dense pair plasma in anti-GJ flow with j m = — 1.5/gj- The same quantities are plotted as in Fig. 14. 


these protons never exceeds ~0.2«gj> at least 10 times less than the 
number density of electrons and positrons; protons do not play a 
substantial role in the electrodynamics of discharges. The number 
density of protons pulled from the NS while electrons and positrons 
are still present near the NS surface depends on the numerical res- 
olution and, hence, our simulations are inconclusive about whether 
this effect is real. However, there is a stage in cascade development 
when the region near the NS gets depleted of pair plasma; at those 
times proton extraction from the NS is certainly a real effect. The 
number of particles left from the previous burst of pair formation 
near the NS is not enough to support the imposed current density 
and screen the electric field until the new portion of the pair plasma 
arrives. A gap with the electric field appears at the NS surface - see 
snapshots at t = 6.730, 6.830 in Fig. 16. This gap is best visible on 


phase portraits of positrons as a second region (the one close to the 
NS) depleted of those particles. The electric field in the gap starts 
accelerating protons from the NS surface, as is clearly visible in 
phase space portrait for protons in Fig. 16. In our simulations this 
second gap does not become large enough to accelerate electrons 
to energies sufficient to start another burst of pair formation near 
the NS surface, but we cannot exclude this possibility in all cases. 
Finally, the main gap reaches the star and pulls out protons as well 
- see snapshot at t = 6.930 in Fig. 16. 

6.2 Flow with j m Ugj > 1 

In this section we describe cascade development for the case when 
the imposed current density is super-GJ. As an example of such 
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Figure 17. Momentum distributions of particles at three moments of time for cascade with anti-GJ current density j m = — 0.5/gj- Positron distributions are 
shown by solid blue lines, electron distributions by dashed red lines, distribution of pair producing photons by dotted black lines. Plots in the top row show 
distributions for particles moving towards the magnetosphere (p > 0) and plots in the bottom row show the distributions for particles moving towards the NS 
( p < 0). Each column corresponds to the same moment of time shown above the plots. Plots in each columns are aligned and share the same values of \p\. The 
following spacial regions were used for plotting the average distribution functions: x e [0, 0.7] for t = 8.350, x e [0, 0.75] for t = 8.510 and x e [0, 0.85] for 
t = 8.710. 


t = 6.560 


t = 6.730 


t = 6.930 


dn 





Figure 18. Momentum distributions of particles at three moments of time for cascade with anti-GJ current density j m = — 1.5/gj. Notations are the same as in 
Fig. 17. The following spacial regions were used for plotting the average distribution functions: x e [0, 0.75] for t = 6.560, x e [0, 0.88] for t = 6.730 and x e 
[0, 1] for t = 6.930. 


flow we consider the case with j m = 1.5/gj- To support this current 
density electrons must move up, towards the magnetosphere, and 
positrons down, towards the NS. There is a continuous source of 
electrons at the surface of the NS, x = 0. Positrons appear during 
the bursts of pair formation and there is no source of positrons 
during the quiet phase, between successive bursts of pair forma- 


tion. To keep the electric field screened there must be both elec- 
trons and positrons present at each point - electrons with larger- 
than-GJ number density moving up and positrons moving down. 
Positrons compensate for the larger-than-GJ number density of 
electrons keeping the total charge density equal to the GJ charge 
density. 
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Figure 19. Screening of the electric field and formation of superluminal electrostatic wave for cascade with / m = — 0.5/gj- There are six snapshots for the 
electric field E, the total change density rj and the charge density of electrons (negative values, blue line) and positrons (positive values, red line) r)±. All 
quantities are plotted as functions of distance x for the part of the calculation domain with intense pair formation. Snapshots are taken at equally separated time 
intervals. Plots in each column are aligned and share the same values of x. The same normalizations for physical quantities are used as in Fig. 1 1 . The three 
thin red vertical lines in each plot mark fiducial points moving with the speed of light towards the magnetosphere. 


Positrons are leaving the domain and as there is no source of 
positrons, some time after the burst of pair formation the domain 
becomes depleted of positrons. This depletion starts at the upper 
end of the domain - at the ‘magnetosphere’ - as positrons are 
moving towards the NS. In the region depleted of positrons electrons 
produce current density [/j < [/ m | and the charge density \ r) | > |t)gjI, 
because their number density is set at regions closer to the NS where 
positrons are still present. This gives rise to the accelerating electric 
field which grows as the remaining positrons are moving towards 
the NS (see Fig. 20). For time shots at t = 9.380-9.500 on the plots 
for the charge density r/ and the current density j there are small 
jumps in both rj and j at the point where the number density of 
positrons drops to zero. The electron number density remains the 
same, but positrons are leaving, enlarging the region depleted of 
positive charges. 

As the gap grows, electrons are accelerated up to higher and 
higher energies and start emitting pair producing gamma-rays. In 
this case, electrons extracted from the NS surface, with the number 


density of the order of \j m \/e, are the particles which ignite the next 
burst of pair creation. As electrons are moving up (as do the first 
pair producing photons), the first pairs are produced at the farthest 
distance where pair formation starts to be possible, near x = x B = 
0.7 L in our example. The injected pairs are picked up by the very 
strong electric field and are accelerated to high energies in less time 
than the primary electrons. Soon they start emitting pair producing 
capable gamma-rays which decay into pairs, see snapshots for t > 
9.420 in Fig. 20. Electrons and positrons are moving in opposite di- 
rections; the pair plasma becomes polarized and starts screening the 
electric field, see plots for r] ± , rj and E at t = 9.500, 9.520 in Figs 20 
and 21. 

The screening of the electric field proceeds similarly to that in 
the case of anti-GJ flow with j m = — 1.5/gj described in Section 6.1. 
The electric field is being screened first near x = x B ; this creates 
a finite size gap with accelerating electric field which moves to- 
wards the NS. The bulk motion of positrons left from the previous 
burst of pair formation is relativistic and the newly created pairs are 
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Figure 20 . Ignition of pair formation in super-GJ flow with j m — 1 . 5 / gj . The same quantities are plotted as in Fig. 11 . 


relativistic too, so the boundaries of the regions where the electric 
field is screened are moving in the same direction (towards the NS) 
with the same speed. The region in between with a non-screened 
field - the gap - is therefore also moving, retaining its size. Elec- 
trons are continuously extracted from the NS to sustain the imposed 
current density. They enter the gap, are accelerated and emit pair 
producing gamma-rays. Pairs are produced by (i) primary electrons 
extracted from the NS surface accelerated by the moving gap and 
(ii) secondary positrons moving towards the NS accelerated in the 
initial screening event. The beam of primary electrons accelerated in 
the gap is clearly seen on plots for particle momentum distribution 
(Fig. 23) as a spike in the electron distribution function for positive 
values of p (the upper plots in Fig. 23). In contrast to the case with 
j m = — 1 . 5 / gi , where the number of particles which were continu- 
ously accelerated in the gap was small and most of the pair were 
produced with momenta directed towards the NS, here more pairs 
are injected with momenta directed towards the magnetosphere - 
see Fig. 23 and plots for particle number density tj± at t > 9.640 in 


Figs 21 and 22. As in the case of anti-GJ flow secondary particles 
have broad momentum spectra - the low energy plasma components 
are present at all stages of cascade development (Fig. 23). 

The pair creation ends when the down flow component of the 
pairs reaches the stellar surface - the gap closes and the accelerat- 
ing electric field is poisoned throughout the strong magnetic field 
region, while the beam component that emits the primary pair creat- 
ing photons and its daughter pairs move up at the speed of light, into 
the magnetosphere, as can be seen in the e + and e~ phase portraits 
with advancing time in Figs 21 and 22. Positrons’ bulk motion is 
towards the NS and some time after the plasma burst leaves - the 
wake of the departing burst being the source of positrons - the upper 
regions become deplete of positrons and the cycle starts again. 

Both in the case of super-GJ flow with constant GJ charge den- 
sity described here and in anti-GJ flow described in Section 6.1, 
the cascade starts at the upper boundary of the strong field region 
when the upper regions become depleted of positrons or electrons 
correspondingly, as the latter move down. The rate at which these 
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Figure 21. Screening of the electric field in super-GJ flow w 


thf m = 1.5/gj. The same quantities are plotted as in Fig. 11. 


regions become charge depleted depends on the imposed current 
density; the higher the absolute value of j m the faster is the bulk 
motion of the charge component moving towards the NS. For [/ m | 
comparable or larger than L/gj | this bulk velocity is close to c. Thus 
the limit cycle period is larger than the relativistic fly-by time over 
the length of the strong field region. In our simulations, that region 
has size 0.7L = 168 m, with fly by time 0.56 ps. However, in a 
real pulsar, the optical depth for pair creation by curvature photons 
emitted by the beam exceeds unity all the way out to heights com- 
parable to the stellar radius. Then the fly-by and recurrence times 
might be as long ~30 ps, with the intermittent plasma starved re- 
gions taking the form of long filaments. Assessing these structures 
and accounting for competition between neighbouring filamentary 
gaps are intrinsically 2D issues, and will be treated elsewhere. 

Because the cascade starts far from the NS, at the largest dis- 
tance where pair formation becomes possible, in our simulations 
the accelerating electric field is unscreened between bursts of pair 
formation above the pair formation zone. In this regard in the case 


of constant GJ charge density super-GJ flow is similar to the anti- 
GJ flow. Namely another cascade zone in the outer magnetosphere 
might coexist with the polar cap accelerating zone. However, as we 
will show in the next section variation of the GJ charge density with 
the distance can significantly change the behaviour of cascades in 
the super-GJ flow. 

6.3 Effects of spatially varying Goldreich-Julian 
charge density 

In this section we address the effect of the mismatch between the 
charge density of the space charge limited beam component and 
the local GJ charge density on the cascade development. In the 
widely adopted steady flow model of polar cap cascades (Arons & 
Scharlemann 1979; Muslimov & Tsygan 1992) this mismatch is the 
reason for the appearance of the accelerating electric field. 

We performed simulations with different scaling of the GJ 
charge density, both in shape (linear, power law, exponential. 
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Figure 22. Filling of computation domain with dense pair plasma in super-GJ flow with j m = 1.5/gj- The same quantities are plotted as in Fig. 11. 


decreasing/increasing with the distance) and in amplitude for all 
the cases considered in the previous sections. For the space charge 
limited flow with 0 < j m /joj < 1 and j m / jai < 0 we found no 
qualitative difference in cascade behaviour due to variation of the 
GJ charge density. Namely, for anti-GJ flows, < 0, the re- 

gions far from the NS became charge-depleted first and the cascade 
started at the largest distance where pair formation is possible, with 
the acceleration zone propagating towards the NS. Sub-GJ flows, 
with 0 < jm/jd < 1, - if the imposed current density was less 
than the local value of the GJ current density r]Q S (x)c everywhere 
- remained low energetic and had two components: a moderately 
relativistic beam propagating through a cloud of trapped particles. 
The main properties of space charge limited flow for these cases 
described in Sections 5 and 6.1 are not affected by variations of the 
GJ charge density with the distance. 

The flow with super-GJ current density j m /joi > 1, however, can 
be strongly affected by variations of the GJ charge density - we 
observed different flow behaviour when the absolute value of the 


flai/B increased with the distance from the NS. In that case the gap 
can appear close to the NS surface, in contrast to all other cases 
when it appeared at large distances. If such gaps can generate pairs, 
the regions above such gaps, far from the NS, remain filled with 
dense pair plasma at all times. 

Below we describe in detail an example of such flow using simu- 
lations with exaggerated charge density contrast in order to clearly 
demonstrate the above mentioned effects. We assume that the GJ 
charge density changes with the distance x as 

n{x) = riai [l +a (0] , (16) 

where r/Qj is the GJ charge density at the surface of the NS and a is 
a positive number. 

Very close to the surface (x <5C R*) all the sources of inhomo- 
geneity of >Jgj/ 5 can be approximated in this manner. For example, 
if variations of the GJ charge density are because of the field line 
curvature (Arons & Scharlemann 1979) the parameter a in equation 
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Figure 23. Momentum distributions of particles at three moments of time for cascade with super-GJ current density j m = 1.5/gj- Notations are the same as in 
Fig. 17. The following spatial regions were used for plotting the average distribution functions: x e [0, 0.795] for t = 9.640, x e [0, 1] for t = 9.880 and x e [0, 
1] for t = 10.100. 


(16) for the dipolar magnetic field would be given by the following 
expression 

3 0, cos />* sin x L 

uas ^ , (17) 

4 cos / — (3/2) 0„ cos </>* sin x R* 

where 0* and are colatitude and azimuth of the field line at the 
NS surface [cf. equation (12) in Arons & Scharlemann (1979)]. 
Note that for favourably curved magnetic field lines cos 0* < 0 and 
so a as > 0. If GJ charge density variations are due to inertial frame 
dragging (Muslimov & Tsygan 1992) then 

L 

flMT 3/C g cos x — , (18) 

K g = 2 GI / R^c 2 ~ 0.15 and / = NS’s moment of inertia [cf. equa- 
tion (32) in Muslimov & Tsygan (1992)]. 

Physical parameters of our simulations are chosen in such a way 
that the length of the accelerating gap is smaller than that of xb- 
This setup differs from those described in Sections 6. 1 and 6.2 in 
that the vacuum potential drop over the domain is larger — AV = 
9 x 10 14 V. 8 The parameter in expression (16) for the GJ charge 
density is a = 0.7, so that the GJ charge density varies linearly from 
— |^qj| to — 1 -7| » 7 qj | throughout the domain. 9 The imposed current 
density is j m = 2 /<?■ j = 2 ))qjC, so that everywhere in the domain 
the imposed current density is larger than the local value of the GJ 
current density, j m > rjaj(x) c. As the variation of the GJ charge 
density with the distance (after accounting for decrease in S; see 
the second paragraph of Section 6) is not larger than 15 per cent 
(Hibschman & Arons 2001), our setup can be considered as a toy 
model for cascades in young pulsars with j m > 1.15 /gj. All other 
parameters of the model are the same as those of the models in 
Sections 6.1 and 6.2. 

8 See Footnote 5. 

9 This amount of variation is huge compared to what actually occurs - 
for example, if the variation is due to inertial frame dragging, for realistic 
parameters a ~ 0.01, if L ~ r pc . 


In contrast to all other cases considered in previous sections the 
dense pair plasma is always present at large distances from the NS - 
the intermittent temporal gaps separating periods when the domain 
is full of plasma disappear. So, we illustrate our example with two 
series of snapshots in Figs 24 and 25. 

When the whole domain is filled with pair plasma the current 
density is constant, / = j m . The number density of electrons is larger 
than [ tjGj|/e and positrons (which are left from the previous burst of 
pair formation) are necessary to keep the charge density equal to the 
local value of the GJ charge density. The absolute value of the GJ 
charge density is smaller closer to the NS surface, rj gj is negative 
and more positrons are necessary to screen the electric field there 
than at larger altitudes. Positrons, on average, move towards the NS, 
and, as they slam on the NS surface and leave the domain, positrons 
from higher altitudes should come closer to the NS in order to keep 
the electric field screened. At some altitude the required flux of 
positrons towards the NS becomes larger than the positron flux that 
can be extracted from an adjacent higher altitude region without 
making that region charge starved - the resulting charge starvation 
causes a starvation electric field to appear close to the NS (see 
snapshot at t = 9.528 in Fig. 24, the gap is best visible in the phase 
portraits on electrons and positrons). In the region above the gap, 
there are still enough positrons to screen the electric field - less 
positrons are needed here as the absolute value of the GJ charge 
density is higher there and so the mismatch in charge density \j m /c 
— r](x) | is smaller. At lower altitudes, after all positrons which had 
been there at an earlier time, before the gap appeared, leave the 
domain, the influx of positrons from above is unable to keep the 
electric field screened and so the gap extends up to the NS surface 
(snapshot at f = 9.636 in Fig. 24). 

The upper end of the gap extends towards higher altitudes as 
positrons move towards the NS, depleting higher altitude regions 
of positive charge. While the gap develops, until the ignition of 
the new cascade, the current density remains close to the imposed 
current density - electrons are being extracted from the NS surface 
and provide the required current density (see plots for j in Fig. 24). 
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Figure 24 . Ignition of pair formation in the flow with linearly varying GJ charge density and the imposed current density j m = 2/qj. The same quantities are 
plotted as in Fig. 1 1 . 


The current density in the domain is super-GJ, j > GJj, and as there 
are too few positrons in the gap the charge density there becomes 
super-GJ, |t?| > I ?7 gj I i this gives rise to a strong accelerating electric 
field in the gap (see the plot for E at t = 9.636 in Fig. 24). Both 
electrons extracted from the NS surface and positrons flowing from 
higher altitudes are being accelerated in the gap and as the gap 
widens their maximum energies get higher and they become the 
primary particles igniting the new discharge cycle (snapshot at t = 
9.744 in Fig. 24). 

Filling of the gap with pair plasma proceeds quickly, as the pairs 
are injected throughout the whole length of the gap (time shots at 
t = 9.780and9.816 in Fig. 25). This happens because both elec- 
trons and positrons, moving in opposite directions, are emitting pair 
creation capable photons, and their number densities are compara- 
ble. The gap is filled with plasma before the particles created in the 
previous burst of pair formation leave the higher altitude region, 
and so, in contrast to the cases described in Sections 6. 1 and 6.2 the 
higher altitude regions are always filled with dense pair plasma. 


Screening of the accelerating electric field during each burst of 
pair formation proceeds slightly differently than in cases described 
in the previous sections. Particles are injected throughout the whole 
length of the gap and not at the gap’s one end. Simulations shown in 
Fig. 25 have two ‘centres’ within the gap where screening starts. In 
Fig. 26 we show in more detail than in Fig. 25 how the accelerating 
electric field is being screened by newly created pair plasma about 
the time t = 9.816. Fluctuations of charge and particle number 
densities are significantly less pronounced than those in Fig. 19 
because the accelerating field was due to mismatch of the charge 
density \j m /c — t?| < | G J ? 7 1 . The region of the screened electric field 
spreads from two centres towards the edges of the gap. In Fig. 26 
we show spreading of the second zone with the centre around x — 
0.2; the first ‘centre’ was at x — 0.1 where the screening started 
earlier and the electric field is now almost completely shorted out. 
This spreading occurs again in the form of an electrostatic wave 
phase velocity of which is larger than c and amplitude decreases 
with time. 
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Figure 25. Screening of the electric field and final stages of pair formation in the flow with linearly varying GJ charge density and the imposed current density 
j m = 2/gj- The same quantities are plotted as in Fig. 11. Note that the time interval between the first two snapshots is smaller than those between the rest of 
the plots. 


The position where the gap appears depends on the imposed 
current density and on the variation of the GJ charge density - for 
higher/ m and/or stronger variation of |GJi;| the gap starts developing 
closer to the NS. The accelerating electric field in this gap is smaller 
than that in case with the constant GJ charge density from Section 
6.2 because the number density of positrons is high and so not the 
whole charge density j m /c is contributing to the electric field. If this 
electric field is not enough to produce pairs, then the gap grows 
up to high altitude until most positrons leave the domain; then a 
vacuum-like gap, with the charge density of the order of j m /c, starts 
developing at the outer end of the domain as in the case of constant 
GJ charge density from Section 6.2. 

The reason why variation of GJ charge density does not affect 
space charge limited flow with sub-GJ current density is the pres- 
ence of the trapped particle ‘cloud’ component which can adjust to 
any charge density variation (assuming the current density remains 
sub-GJ), as was mentioned in Section 5. For flow with anti-GJ 


current density the accelerating electric field appears in the region 
depleted of electrons, which on average move towards the NS. In 
this case regions closer to the NS always have the larger number 
density of particles of the same sign as tjgj and will become charge 
starved last, and so the gap develops at higher altitudes. For flows 
with the super-GJ current density when the absolute value of the 
GJ charge density decreases with altitude, a < 0, fewer positrons 
are necessary at lower altitudes to poison the electric field than at 
higher altitudes. Positrons are moving towards the NS and the low- 
altitude region becomes charge starved after all the others and the 
gap develops at higher altitudes as in the case of the flow with a 
constant GJ charge density. 

Therefore, the regions with super-GJ flow can have two qualita- 
tively different behaviours depending on the value of the imposed 
current density and the variation of the GJ charge density. The cas- 
cade either starts close to the NS and high altitude regions are always 
filled with dense pair plasma, or the cascade starts at high altitude 
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Figure 26. Screening of the electric field and the formation of superluminal electrostatic wave in the flow with linearly varying GJ charge density and the 
imposed current density j m = 2/gj. The same quantities are plotted as in Fig. 19. 


and between the bursts of pair formation the high altitude regions 
are charge starved with vacuum-like accelerating electric field - 
in the second case an additional cascade zone might appear in the 
outer magnetosphere. In the first case the period between succes- 
sive bursts of pair formation should be of the order the gap’s fly-by 
time and can be as short as ~r pc /c ~ fractions of microseconds; in 
the second case the repetition rate will be larger than the fly-by time 
of the strong field zone, ~ tens of microseconds. From our 

toy model we cannot make quantitative prediction about the values 
of the imposed current density and magnitudes of the GJ charge 
density variations which separate those behaviours. 

6.4 On stationary cascades 

In all models previously considered in the literature, the size of the 
accelerating gap was limited by the opacity of the magnetosphere to 
the gamma-rays: pairs were injected when the optical depth to pair 
creation reached unity and those pairs screened the electric field. 
In our simulations with the constant GJ charge density, the electric 
field in the gap monotonically increases with distance, and the gap 
grows with time. That electric field could be screened only by the 


injection of particles with the charge sign opposite to the charge 
sign of primary particles, and those particles, when injected, were 
accelerated in the direction opposite to the direction of motion of the 
primary particles by a strong electric field. They start emitting pair 
producing gamma-rays after travelling a short distance, injecting 
pairs and so screening the electric field in their way. That led to the 
motion of the gap as a whole or to the gap closure, if the other end 
of the gap moved with subrelativistic velocity. 

In the case of a flow with varying GJ charge density, described 
in Section 6.3, the initial field was not a monotonic function of the 
distance, but the gap was growing and, again, the only way to limit 
the gap and screen the electric field was by the injection of particles 
of the opposite sign. Once injected those particles ultimately de- 
stroyed the gap, since their flux exceeded that of the primary beam 
- stationary equilibria with counter streaming particles of the op- 
posite sign exist only when the counter streaming flux is less than 
the primary flux, as in the Arons & Scharlemann (1979) model. 
Stationary particle acceleration and pair creation would be possible 
if pairs are injected (mostly) outside the gap, so that there will not 
be too many particles with charge sign opposite to that of primary 
particles within the gap. 
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For example, the stationary flow model of Arons & Scharlemann 
(1979) constructs a finite length acceleration zone bounded by the 
stellar surface below and the abrupt onset of pair creation above, 
which screens the electric field in a thin layer [the pair formation 
front (PFF)]. The space charge limited (electron) beam extracted 
from the surface provides the primary particles that radiate pair cre- 
ating gamma-rays. Almost all of the pair creation occurs above the 
PFF. A small amount of trapping within the transitional PFF layer 
sends charges of sign opposite to the primary beam (positrons) back 
down towards the surface with flux small in comparison to that of 
the primary beam, while extra electrons are added to the beam, 
restoring the charge density to equality with the Goldreich-Julian 
density. This restoration causes the screening (‘poisoning’) of the 
accelerating field. In that model, E — > 0 at the PFF is imposed as 
a boundary condition, using the conclusion that dense pairs ( n± 
lhGj/e|) flash into existence at a rather well-defined height - a 
physically correct conclusion when curvature radiation dominates 
the gamma-ray emission and one photon absorption in the magnetic 
field is the major opacity source, since each primary beam particle 
emits many pair creating gamma-rays and opacity is a rapidly vary- 
ing exponential function of the photon energy. The application of 
that boundary condition, along with the space charge limited flow 
condition at the NS, has the consequence that the charge density 
of the beam has a unique value tjbeam- Since j = crj b eam i this model 
becomes discrepant with the magnetospheric current density j m , in 
general. 

We tried to find an Arons-Scharlemann like solution with the 
acceleration gap limited by the PFF. For the linearly varying GJ 
charge density given by equation (16) and model parameters from 
Section 6.3 we explored the parameter space by running differ- 
ent simulations for different values of the imposed current density 
starting from j m = 2 j and decreasing it up to j m ~ y'gj when pair 
production ceased. In each subsequent simulation with smaller j m 
the position of the first pair injection moved further from the NS 
because the accelerating electric field was also smaller. In all cases 
when pair injection occurred inside the gap the gap boundary starts 
moving - we never saw a stationary PFF (nor did we see a stationary 
PFF in any other simulations). So, the conclusion implied by the 
simulations is that the time asymptotic state is not a steady flow 
similar to Arons & Scharlemann’s. Even though acceleration in the 
charge starved gaps is limited by pair creation with thin boundaries 
between pairs and quasi-vacuum, reminiscent of Arons & Scharle- 
mann’s PFFs, the gaps move, either up into the magnetosphere 
or down towards the star - there is no truly stationary state, and 
fully developed limit cycles appear to be the answer, at least in one 
dimension. 

We now exhibit a novel steady flow model, which takes advantage 
of the properties of non-neutral beam flow when rjaj/B is non- 
uniform, that does not require poisoning by the pairs to produce 
E = 0 at the gap’s top. This model exploits the variation of the 
charge density with distance, in the case when the current density 
is larger than the local GJ current density only up to some height 
h s , then at distances larger than h s - where the non-neutral flow 
becomes sub-GJ - the electric field will be screened on the scale of 
the local Debye length, as discussed in Section 5. In such a situation 
the accelerating electric field exists only up to the distance h s , above 
which it is shorted out by a mixture of trapped electrons - those 
with two turning points in their orbits - and free pairs. 

We assume the GJ charge density varies linearly with the distance, 
equation (16), as discussed in Section 6.3. By hypothesis, the flow 
is steady and the current density is equal to the imposed current 
density j m = fy' ( ? J , f > 1 (see equation 12). If f — 1 < a, at 


some point the non-neutral beam flow becomes sub-GJ. From the 
Poisson equation for the electric field (6), in the absence of pairs 
we have 

E s = 47Tt?Gj£ 

At the distance 
§ - 1 

h s = 2- L, (20) 

a 

E s changes sign as the flow becomes sub-GJ. Identifying h s with 
the gap height yields the gap’s potential drop to be 

8 n 9 (f - l) 3 

AV s , gap =-7i^ J L 2 ^-^, (21) 

increasing in proportion to the imposed current density. If the re- 
sulting potential drop in the gap is large enough for particles to 
inject pairs within the gap with a number large enough to form a 
positron back flux larger than the primary electron flux, the flow 
will be non-stationary, like that described in Section 6.3. But if the 
potential drop is such that copious pairs will be injected at distances 
larger than h s , positrons will be subjected only to a fluctuating, rel- 
atively small electric field, which sustain the cloud component at 
altitudes higher than h s , where the flow is sub-GJ. Only a fraction 
of the positrons will be adverted into the gap, and, if this fraction 
will be much smaller than the GJ charge density the flow in the gap 
will be not disturbed enough to make it non-stationary. 

As an example of a flow with E(h s ) = 0 in the non-neutral current 
flow, pair creation and non-neutral trapping, we show snapshots 
in Fig. 27 of the flow with the same parameters as in Section 6.3, 
except with a carefully chosen current density which was set large 
enough to allow pair formation within the domain (at x < xb), but 
small enough so that the particle back flux does not destroy the gap. 
The imposed current density in this model is j m = 1 .059 j. Above 
the altitude x = h s 0. 1 7L the flow is sub-GJ and clearly has beam- 

cloud structure qualitatively similar to that in sub-GJ flows from 
Section 5 (see phase portraits for electrons). The charged cloud 
of trapped electrons screens the electric field at x > h s . Primary 
electrons, when accelerated up to high energies, emit gamma-rays 
which decay into electron-positron pairs above the gap. 

Both secondary electrons and positrons mix with the cloud com- 
ponent and some positrons are adverted towards the gap. The gap 
shrinks a bit and the resulting potential becomes smaller than a 
critical value necessary to sustain pair production. When the num- 
ber density of positrons drops, the gap gets bigger and the pair 
formation starts again. The gap, however, never disappears and the 
cascade is close to a stationary configuration. This configuration is 
quite sensitive to the imposed current density; large gap fluctuations 
appear at imposed current densities only a few per cent larger than 
j m = 1.059_/gj. Our model, however, has an exaggerated charge 
density gradient and a very high voltage, and we expect that the 
mixing of positrons and their advection for real pulsar parameters 
would be less efficient and, hence, there might be a larger parameter 
range (i.e. an interval of the imposed current densities) where this 
gap plasma flow can sustain quasi-stationary cascades. 

Such stationary configurations require the imposed current den- 
sity to be larger that the local GJ current density at the NS surface 
but then becoming less than the local GJ current density at some 
distance from the NS. Variations of the GJ charge density (after 
accounting for the decreasing magnetic field) would not be larger 
than ~15 per cent (Muslimov & Tsygan 1992; Hibschman & Arons 
2001) and, hence, the imposed current densities for stationary cas- 
cades are within the range |_/ m — y'^jl < 0.15, jm/ jos > 1. Taking 
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Figure 27. Plasma flow with quasi-stationary cascade. The imposed current density is j m = 1 ,059/gj; it is ‘fine-tuned' in such way that the accelerating electric 
gap near the NS surface does not disappear. All other physical parameters are the same as in the case shown in Figs 20 and 2 1 . The same quantities are plotted 
as in Fig. 11. Note that the time intervals between snapshots are not equal. 

into account that the resulting gap should not be very large to prevent 
pair injection within the gap, that range is even smaller. 

Fig. 1 makes clear that in the force-free theory, the current density 
does not mostly fall in this narrow range, lending support to the 
conclusion that limit cycle and trapping behaviour are the generic 
physics for polar flow, with and without pairs. 

7 DISCUSSION 

These calculations, and those reported in Timokhin (2010), were 
designed to investigate the theoretical issue of how the pulsar mag- 
netosphere with current flow determined by the magnetosphere’s 
global dynamics couples to the NS through the polar cap beam 
acceleration and pair creation zone. This problem was investigated 
30-40 years ago, primarily using the order of magnitude estimates 
and analytical, steady flow models of charge-separated beams; the 


accelerator was modelled as having voltage fixed by pair creation 
and geometry. That led to a determination of the charge density 
in the accelerator with the current then simply being charge den- 
sity times the speed of light. 10 No effort was made to show that 
the current density estimated actually matched that required by the 
global structure - the models yielded currents with the correct or- 
der of magnitude, and in the absence of actual solutions for the 
global structure of the magnetosphere, that answer was regarded as 
good enough, although scepticism was expressed that operating the 
accelerator model as having fixed voltage could actually produce 
the correct answer (see Mestel 1999 and references therein). The 
calculations reported here show that indeed when the accelerator is 


10 In differing degrees, the charge density was made up of counter streaming 
relativistic beams, leading to current density lying between t]gjc and 2r/ajc. 
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operated with current fixed rather than voltage, the polar cap accel- 
erator does work, and works in a fully time-dependent manner, as 
first conjectured by Goldreich & Julian (1969) and Sturrock ( 1971), 
but with behaviour different from previously published models. 11 

Our results show that when the acceleration and cascades are one 
dimensional, acceleration and pair creation exhibit limit cycle be- 
haviour, with cycle time somewhat larger than the relativistic fly-by 
time over the length of the accelerator. When the ID approxima- 
tion is realistic (very young pulsars), the length of the accelerator 
is less than the polar cap diameter r pc & 144P -1 ^ 2 m, suggesting 
quasi-periodicity in the acceleration and pair creation on time-scales 
> 10 -6 P -1 ^ 2 s. As was suggested long ago, such variations, if re- 
flected in the radio emission, might be the origin of the microstruc- 
ture in the radio emission (Ruderman & Sutherland 1975). If limit 
cycle behaviour is robust with respect to multi-dimensional con- 
siderations, this dynamical behaviour might turn out to be a useful 
model for microstructure (in contrast to non-linear optical phenom- 
ena in the transfer of radio waves through the pair plasma, for 
example). 

When charges are fully bound to the surface (solid surface with 
high binding energy and no atmosphere; Medin & Lai 2010), cur- 
rent flow always adjusts to the magnetospheric load through pair 
creation discharges, as was shown in Timokhin (2010). The lowest 
energy particles in the pair discharge exhibit transient trapping in the 
fluctuating electric field. Those reversals of particle momenta allow 
both current and charge densities to adjust to the presumed force- 
free conditions, for any value of the magnetospherically imposed 
current density. 

At long periods/weak B, pair creation becomes impossible, 
the magnetosphere is a vacuum with no conduction current flow 
and spin down occurs through the generation of strong vacuum 
waves (Pacini 1967). Pulsars enter this regime by crossing the 
pair creation ‘death valley’, theoretically expected for magne- 
tospheric voltages <J> m = ,/Wr/c ~ 10 13 (P/IO -15 ) 1 / 2 /* -3 / 2 V < 
10 12 = $ death V [a restatement of the death line based on cur- 
vature gamma-ray emission and one photon magnetic conversion 
estimated by Sturrock (1971)]. The fact that this death line describes 
the disappearance boundary of radio pulsars in the P, P diagram 
reasonably well (see e.g. fig. 15.1 in Arons 2009) provides a strong 
hint that low-altitude, polar cap pair creation with high-energy beam 
acceleration and one-photon magnetic conversion of the curvature 
y-rays has something to do with the pulsar activity. 

A warm, quasi-neutral atmosphere plasma atmosphere overlying 
the magnetized ocean and solid crust, with no restrictions on charged 
particle motion along B, is the likely surface state of most pulsars, 
either because of the residual heat of the NS or because of polar cap 
heating by precipitating particles from the magnetosphere (from 
local pair discharges or from the return current flow). The upper 
atmosphere can then freely supply charge, in a manner very similar 
to a space charge limited current flow from the cathode of a classical 
vacuum tube. Our simulation results on space charge limited flow 
with (and without) pair creation reveal a variety of noteworthy new 
features. 

We have made the first determination of the polar accelerator’s 
behaviour under conditions where the current density, not the volt- 


11 Following Shibata (1991), when the pulsar is not near ‘death valley’ 
in P — P space, the load inserted into the magnetospheric circuit by the 
pair creation discharges creates a small perturbation of the magnetospheric 
current j m , thus justifying our treatment of j m as an external parameter in 
the description of the discharges. 


age, is held fixed. We found that space charge limited flow is not 
always high-energy. Emission of pair creating gamma-rays does 
not occur on all polar field lines, even when <t> m > $ deat h. The 
force-free model of the magnetosphere suggests the polar current 
flow includes sub-GJ flow regions with 0 < j/jai < 1, where the 
low energy current flows (with the Lorentz factor /beam 3) with 
no (curvature) gamma-ray emission and no pair creation (inverse 
Compton gamma-ray emission and y—y conversion to pairs are both 
negligible), with adjustment of the current density and charge den- 
sity to the force-free values through formation of a trapped particle, 
non-neutral hanging charge cloud. This kind of quasi-stationary 
flow occupies most of the polar flux tube for the aligned and al- 
most aligned rotator, but progressively disappears as the obliquity 
increases. The force-free model also requires regions of distributed 
return current flow,y / j gj < 0, and, as the obliquity increases, regions 
of super-GJ current flow,_///GJ > 1 > occupying most of the polar flux 
tube as the obliquity approaches 90°. Both the return current and 
super-GJ regimes exhibit unsteady, high-energy current flow (an 
unsteady beam), driving discharges copiously creating pairs whose 
character is similar to that encountered when the atmosphere is ab- 
sent and the surface is a strongly bound solid. These regions, as 
projected on to the polar cap, are shown in Fig. 1. In the special 
case of j m / y'(?j being slightly larger than unity, the plasma flow can 
sustain quasi-stationary cascades, and a new class of such models 
was described in Section 6.4. But generally speaking, such station- 
ary regimes represent a singular case, rarely if ever achievable by 
magnetospheres described by the force-free model. 

The only previous quantitative study of time-dependent cascades 
in space charge limited flow is that of Levinson et al. (2005). They 
used a two-fluid model and concluded that chaotic pair formation 
takes place throughout the whole zone where pair formation is pos- 
sible. Our simulations do not support their conclusions. Particle 
trapping plays a significant role in adjusting of the plasma flow to 
the required current density and, hence, the plasma cannot be ade- 
quately represented as two fluids (electrons and positrons) each with 
its own unidirectional velocity. These two-fluid representations in- 
troduced additional rigidity into the system and this, in our opinion, 
led to the formation of a strong chaotic electric field everywhere in 
the domain. 

The most recent analytical treatment of the problem was pre- 
sented in Beloborodov (2008). Our results support Beloborodov’s 
(2008) general conjecture about the character of plasma flow in 
the force-free regime, that the sub-GJ flow is low-energetic and 
pair formation is possible only in super-GJ or anti-GJ current flow. 
The simulations differ from his qualitative picture of what would 
happen in several respects: (a) the low energy flow with trapped 
particles does not retain the spatial oscillations of the cold flow, 
replacing those by the two-component beam/cloud structure - the 
warm beam has non-oscillatory velocity, while the trapped parti- 
cles, not included in his picture, are those with two turning points 
in their orbits; (b) bursts of pair cascades repeat after longer time 
than h/c , h being the gap height; and (c) in most cases the gap is 
not destroyed but moves as a whole, usually relativistically. 

The discharges were given a novel treatment. After the sem- 
inal efforts of Ruderman & Sutherland (1975) but prior to our 
work, models universally incorporated their assumption that when 
pair creation is copious (many convertible gamma-rays per primary 
high-energy particle, as is the case in the curvature radiation cas- 
cades typical of young, high voltage pulsars), at and above the 
height where pairs start appearing, £][ drops to zero and stays that 
way at all greater altitudes. Arons & Scharlemann (1979) formal- 
ized this into the dynamics of a PFF, showing that when pair creation 
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is copious, the transition between a charged starved region where 
E\\ ^ 0 and the pair-dominated region where E\\ can be set equal 
to zero is thin - that structure necessarily involves one sided parti- 
cle trapping, under pulsar conditions, thus causing the formation of 
counter streaming components in the current flow. Our solutions ex- 
hibit this transition between the dense plasma and the quasi-vacuum 
regions outside. But, since E\\ was obtained from a dynamical field 
equation, we nowhere assumed E \ | had to be zero everywhere out- 
side (above or below) a pair discharge. Instead, the pair creation 
itself creates the polarizable plasma that dynamically resets to 
zero inside the discharge, while outside the electric field was al- 
lowed to float to whatever is dynamically required. The resulting 
models thus exhibit macroscopic intermittency - a pair discharge 
shorts out £|| , and then as the cloud of pair plasma moves out into 
the magnetosphere or towards the surface, a gap reforms in which 
residual particles accelerate, radiate gamma-rays and form a new 
discharge. The fact that pair creation opacity declines with distance 
from the NS was modelled by setting the magnetic field equal to 
zero in the upper part (typically the upper 30 per cent) of the sim- 
ulation domain. That Ansatz forces the one photon pair creation 
opacity to be zero in the upper part of the domain. 

It must be said, however, that the resulting model of the in- 
termittency is very simplified - probably oversimplified. The rep- 
resentation of the ‘magnetosphere’ as a region of negligible pair 
creation optical depth to gamma-rays emitted from just above the 
polar cap is the right general idea, but the ID character of our simu- 
lations almost certainly leads to an overemphasis on the coherence 
of the intermittent discharges - they form coherent slabs of pair 
plasma, separated by coherent slabs of quasi-vacuum. In reality, the 
low opacity region where discharges must end appears at heights 
r pc = polar cap and polar flux tube transverse radius. Then, 
the acceleration zone left behind a discharge created plasma cloud 
is long and skinny, with the narrow ‘wave-guide’ geometry (and the 
physics of the polar flux tube’s boundary) playing an essential role 
in the nature of the quasi-vacuum field. Under these circumstances, 
it is unclear and unknown whether the discharges would preserve 
anything like the coherence exhibited in the results reported here - 
the formation of multiple ‘lightning bolts’ is perhaps more likely, 
simultaneously coexisting within the polar flux tube, the electric 
fields of the discharges affecting the dynamics of their neighbours, 
causing the discharges to influence each other - a complex dynam- 
ical system quite likely to be chaotic. 

In addition, the extent of the gaps formed between discharges 
requires consideration of all the sources of pair conversion opacity, 
not only one photon conversion in the B field. Even when magnetic 
conversion drops to zero, y—y conversion with gamma-rays collid- 
ing with soft photons from the atmosphere (either heated polar cap 
or overall warm surface, if the star is young enough) can provide 
discharge initiating pairs within quasi-vacuum gaps, limiting the 
extent of such gaps to less than what occurs when the opacity is set 
equal to zero. These issues require multi-dimensional modelling of 
the accelerator and the photon transfer, as well as extending the ra- 
diation transfer model, improvements required before the possible 
direct consequences of low-altitude discharges for observations can 
be addressed. Those consequences are the pair multiplicity of the 
outflow, the heating of the polar caps by the discharge components 
that move towards the NS, and the possible direct collective emis- 
sion of photons by the time-dependent currents in the discharges. 

The multiplicity ( k± = « pa ir/«Gj) within individual pair creation 
bursts show the traditional results - when a discharge occurs, k± 
~ ybeam"t ± c 2 /£ C urvature ~ 10 3 Ocurvature = energy of pair producing 
curvature photons), somewhat enhanced by the conversion of the 


synchrotron photons in the cascade - thus, multiplicities within in- 
dividual bursts up to ~10 4 are observed. Intermittency reduces the 
overall multiplicity of the outflow - the magnitude of this reduction 
is hard to judge without a multi-dimensional model. The intermit- 
tency reduction is offset by continued pair creation as the plasma 
clouds and their high energy leading edges drift out to heights of the 
order of the stellar radius and more, emphasizing again the need for 
a multi-dimensional treatment. However, the multiplicity is unlikely 
to be as high as is inferred from PWNe (Bucciantini et al. 2011), 
so long as the simple, star centred dipole B field model is retained. 
More general B fields, with smaller radii of curvature of photon 
orbits with respect to B, can lead to enhanced pair creating opacity 
and pair yield - the offset dipole model (Arons 1998; Harding & 
Muslimov 201 1) is the most plausible magnetic anomaly model that 
can be consistent with the dipolar character of the polar flux tube 
revealed by radio observations (Rankin 1990; Kramer et al. 1998). 
Whether sufficiently large multiplicities can be obtained within the 
geometric constraints obtained from the radio polarimetry is under 
investigation. 

Polar cap heating and consequent soft X-ray emission provide an- 
other possible observable that can constrain the model. In common 
with all discharge models, roughly half the energy in each discharge 
is deposited in the crust below the atmosphere and ocean. If intermit- 
tency is neglected, existing observations provide strong challenges 
to all polar cap discharge models, including ours. Intermittency 
reduces the average energy flux. Whether that effect allows this 
polar discharge model to survive confrontation with observations 
(or, better, prove useful in modelling such observations) remains 
to be studied. Simulating the particle back flux from a discharge is 
a particular challenge, requiring resolving the Debye length in the 
pairs, which rapidly decreases as the pair density grows. 

Our results imply that discharges occur on some of the polar 
field lines for all inclinations. These discharges incorporate time- 
dependent, quasi-coherent currents on microsecond and shorter 
time-scales. It has not escaped our attention that such fluctuations 
might be a direct source of radio emission from the low-altitude 
polar flux tube, a region strongly suggested as the site of the radio 
emission by the radio astronomical phenomenology. Although the 
electric fields in our one-dimensional model are wholly electro- 
static, therefore cannot leave the plasma, in a more general multi- 
dimensional setting which has substantial spatial inhomogeneity 
transverse to B, the field components parallel to B have accompany- 
ing components E± perpendicular to B which make the fluctuations 
fully electromagnetic. Therefore, these electrostatic spectra may be 
representative of electromagnetic fluctuations which can leave the 
plasma, and the pulsar. From the point of view of simulations, a 
multi-dimensional, electromagnetic treatment is required in order 
to investigate this possibility. 12 

In the frame of our model we can provide only estimates of the 
energy available for such directly excited waves. If the electrostatic 
oscillations could form an electromagnetic wave which leaves the 
plasma then the energy carried by the wave would be of the order 
of 


W t 


(E 2 ) 

Sn 


( 22 ) 


12 Fawley (1978) used a linear response analysis to study this possibility in 
an analytic model of discharges when strong surface binding suppresses fee 
emission. The work described in this paper is the first to study fully non- 
linear discharges in the free emission regime. Whether this idea for radio 
emission will lead to a useful model remains to be studied. 


Downloaded from http://mnras.oxfordjoumals.org/ by guest on January 9, 2013 


Pair cascades over pulsar polar caps 3 1 



Figure 28. Estimated energy flux W in plasma waves for sub-GJ flows as 
a function of £ = jm/jGJ- W is normalized to P 2 - 


Both kinds of flow, with or without pair formation, cause fluctuating 
electric field. Fluctuating electric field in the sub-GJ flow turns out 
to be too low to provide the energy even for the radio emission. 
The electric field for the sub-GJ flow shown in Fig. 8 is normalized 
to Eq = | ?7gj I GJi in this normalization the estimates for W r is 
given by 

W r ~ 4.8 x 1(T 9 (E 2 )W mi B{ 2 1 P 2 , (23) 


where E is the normalized electric field and JF md is magnetodipolar 
energy losses 


b 2 r 6 q 4 

u * 


4c 3 


(24) 


In Fig. 28, we plot the estimated values of W r = JJi-/B I2 1 P 2 VF m d as 
a function of £ = y m / ioi ■ It is evident that even for the Crab pul- 
sar, known for its very low radio efficiency in terms of spin down 
energy losses ~10~ 8 , the energy in electrostatic oscillations of the 
cold flow cannot power the radio emission. For plasma flows with 
pair formation, screening of the electric field proceeds similarly to 
cascades in the Ruderman-Sutherland model studied in Timokhin 
(2010). We observed the formation of superluminal plasma waves 
during discharges in space charge limited flows, similar to what was 
found for strongly bound surfaces. The range of plasma oscillation 
wavelengths was rather broad, similar to what was observed in dis- 
charges discussed in Timokhin (2010). The electric field for plasma 
flows with pair formation discussed in Section 6 is normalized to 
Eq = \r)Q]\nr vc \ in this normalization the estimates for W r are given 
by 


W T ~ l -(E 2 )W mi . (25) 

The amplitude of the (normalized) oscillating electric field during 
the screening phase of cascade development shown in Figs 19 and 
26 is ~0. 1-0.01 and the energy in such oscillations is more than 
enough to power the radio emission. So, from the energetics point 
of view, it seems that the plasma flows with pair formation studied 
here are candidates for powering the pulsar radio emission through 
direct radiation by the discharge currents. 

Finally, our work may have implications for the formation of the 
particle accelerators in the outer magnetosphere required to account 
for the gamma-ray pulsars. Low voltage, sub-GJ current flow may 
have a Holloway (1973) ‘problem’; in that the non-neutral cloud 
and low energy, non-neutral current-carrying beam cannot cross the 
null surface where ft • B = 0 - in contrast, field lines populated 


with dense plasma from discharges have no difficulty with plasma 
flowing across the null surface, the pair plasma clouds are quasi- 
neutral and easily adjust to the locally required charge density. Field 
lines passing from the low-altitude trapped cloud to the outer mag- 
netosphere that go through the null surface might form a physically 
self-consistent gap reminiscent of the earliest proposals for ‘outer 
gaps’, with the potential for outer magnetosphere discharges that 
could provide a model for the observed gamma-ray emission. This 
issue also requires multi-dimensional modelling. The caustic for- 
mation exhibited by models for gamma-ray pulses suggests that if 
such gaps can form, they are localized (by pair creation?) to regions 
close to the flux bundles where return currents flow. 

Another possibility for particle accelerators in the outer magne- 
tosphere exists in the regions carrying the return current. At some 
time between the bursts of pair formation in the polar cap the elec- 
tric field at large latitudes cannot be screened by pairs created in 
the polar cap. This field can accelerate particles and give rise to 
pair creation via y—y process. As we pointed out at the end of Sec- 
tion 6.1, an outer magnetosphere cascade zone might exist along 
the magnetic field lines carrying the return current. This might be 
a very intriguing possibility in view of the observational evidence 
that the spectrum of pulsar gamma-ray emission does not have a 
super-exponential cut-off, which should be present if there were ab- 
sorption of gamma-rays in the magnetic field. Moreover, modelling 
of pulsar gamma-ray light curves suggests that the gamma-ray emis- 
sion originates from the regions close to the boundary between the 
open and closed magnetic filed lines (e.g. Bai & Spitkovsky 2010). 
The return current regions for a broad range of pulsar inclination 
angles are close to that boundary or/and are within it, e.g. in Fig. 1 
the return current is flowing in and around the current sheet for 
pulsar inclination angles a < 30°. 


8 CONCLUSIONS 

Our principal conclusion is simple - pair creation can occur at 
pulsar polar caps, but (almost) always in the form of fully time- 
dependent current flow (microsecond time-scales for the variabil- 
ity). That time dependence with pair creation allows the current 
to adjust to any magnetospheric load, while simultaneously allow- 
ing the charge density to adjust to the requirements of the force- 
free magnetosphere. We have also shown that a substantial fraction 
of the open field lines (fraction decreasing with increasing obliq- 
uity) solve the current flow problem with a low-energy, non-neutral 
beam carrying the current co-existing with a non-neutral, electri- 
cally trapped particle cloud. This is an essentially time-independent 
local solution. Exploring the consequences of these new results for 
global theory and observations requires extending the calculations 
to multi-dimensional, electromagnetic models for the accelerating 
electric field, and perhaps to background magnetic field models 
more general that the star centred dipole geometry (used here in the 
choice of the magnetic radius of curvature that enters into the pair 
conversion opacity). 
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APPENDIX A: ONE-DIMENSIONAL 
ELECTRODYNAMICS OF THE POLAR CAP 

Let us decompose the magnetic field in the acceleration zone into 
two terms B = Bq + SB, where B 0 is the global time-averaged 
field (the ‘external’ field for our local problem; = B magnetosphere) and 
52? is the fluctuating magnetic field due to local currents in the 
system caused by charge motion in the acceleration zone. 


In our ID approximation the characteristic size of the acceleration 
zone in longitudinal direction, along magnetic field lines, / is much 
smaller than its characteristic width w (the latter being of the order of 
r pc ). The characteristic time-scale of electromagnetic field variations 
due to relativistic charge motion is r ~ l/c, as charges move along 
magnetic field lines. The characteristic scale of the global magnetic 
field variation in longitudinal direction 7 mag is much larger than both 
/ and w. The relation between these scales is / <$( w <5C / mag . 

Electromagnetic field can be determined from Ampere’s and 
Faraday’s laws 

471 1 9 E 

V x B 0 + V x SB = — j H — (Al) 

c c dt 


1 dSB 

V x E = — . (A2) 

c 9 1 

Combining Ampere’s (Al) and Faraday’s (A2) laws by eliminating 
the electric field we get an equation for the magnetic field 8B . 

1 d 2 SB , 4tt 

-V 2 8B=-Vxi + V-B „. (A3) 

c l at- c 

Only the perpendicular components of the magnetic field Bo t ±, 8B± 
affect the accelerating electric field E y . Now using the perpendicular 
component of equation (A3) we estimate 8B±. 

The estimates of each of the terms in equation (A3) are as follows. 
For terms with 8Bi we have 

8B ± 

~P~ 


(y 2 8B) ± r 
and 

1 d 2 8B 
c 2 dt 2 


8B ± 8B ± 

l 2 w 2 


(A4) 


8B ± 

(ct ) 2 


8B_ l 


(A5) 


The perpendicular component of the global magnetic field changes 
in the longitudinal direction on the scale / mag , in the lateral direction 
it changes on the scale w and, hence, 


(V 2 Bq) j 


Bo,± B 0 ± 

/2 yjl 

mag w 


Bo,i_ 

m2 


(A6) 


Current flows along magnetic field lines and its perpendicular com- 
ponent is of the order of (Bo, ±/Bo)j ~ ( 1/ pc)j , and as the radius of 
curvature of magnetic field lines p c is much larger than any of the 
scales in our problem we have 

j 7 1 / l j j 

(V x j) ± ~ — — —r~- i ~ ' (A7) 

W I W p c l W 

Combining estimates (A4)-(A7) from equation (A3) we have 

2 


8B_ l 
w 


471 


- —J + 


B, 


0 , 1 . 


The order-of-magnitude version of Ampere’s law (Al) is 


Bo,± + 8 

w w 


An . En 

— j 3 > 

C CT 


(A8) 


(A9) 


and from equation (A8) it follows that the fluctuating magnetic field 
due to charge motion introduces only a second order term in l/w 
and it can be neglected. So, the Ampere law in one dimension has 
the form 

- - 471 (; - J^V x B 0 )||) - ~An(j - j m ) . (A10) 

The same expression for the accelerating electric field through the 
current density as in equation (A10) can be obtained from the Gauss 
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law and charge conservation [see e.g. Timokhin (2010), Appendix 
A]; in that case j m emerges as an integration constant corresponding 
to the time-average current density flowing trough the system. In 
one dimension the system is essentially electrostatic; charges cre- 
ate only electric field, and naturally both Ampere’s law and Gauss’ 
law reduce to the same equation. The magnetic field perpendicu- 
lar to the background stellar B field is negligible, both due to the 
background magnetospheric current and due to the rapidly variable 
currents within the acceleration zone considered in this study - 
in particular, to the lowest order in ( 1/ p c ) 2 , the particle orbits are 
determined by the stellar field Bo and consideration of the full elec- 
tromagnetism of the discharges is not required to characterize the 
discharge dynamics. 


APPENDIX B: STATIONARY 
ONE-DIMENSIONAL SPACE CHARGE 
LIMITED FLOW 


The first integral of this equation is easily obtained by multiplying 
both part by d y /ds 



- (y - vo) 



(B8) 


Under the conventional assumptions the space charge limited flow 
starts at the surface with zero velocity and the electric field is com- 
pletely screened, so yo = 1 and (dy/d.s)o = 0; with these assump- 
tions equation (B8) becomes 

(£)' = 2 (fU^T- K+1 ). <B9> 

It is more convenient to analyse the properties of the flow in 
terms of the spatial component of the four-velocity (the normalized 
momentum p = yv/c). In terms of p equation (B9) becomes 


B1 General equations 

Let us consider a stationary one-dimensional beam of equal particles 
with mass m and charge q having the same sign as the GJ charge 
density rjoj flowing from the surface of the NS at x = 0 into the 
magnetosphere at x > 0. The current density of the beam is a fraction 
f of the GJ charge density, 

j = Sj tar- (Bl) 


For the energy of particles in the beam we have 

mc 2 y + q<p = mc 2 y 0 + q<po , (B2) 


where y = (1 — ir /c 2 )~ 1|/2 is particle’s Lorentz factor, v is particle’s 
velocity, <p is the electric potential; quantities with the subscript 0 
refer to their values at NS surface, at x = 0. The electric potential 4> 
is given by the Gauss law 


d 2 cp 
dr 2 


= —4n(q - r? G j) , 


(B3) 


where q = j/v is the charge density of the beam. We consider the 
case when v, y, q and <p are functions of distance x, while qaj is 
constant. 

Expressing <p through y from equation (B2) and q through j = 
vq we get an equation for the beam’s Lorentz factor 


d 2 y 

dx 2 


4 nq Gs q 
me 2 


W 2 


(B4) 


The plasma frequency for a mildly relativistic plasma with the GJ 
charge density consisting of particles with charge q and mass m is 


w p,GJ 


( 47t^ G j<? 


1/2 


(B5) 


it) (B10) 

The right-hand site of equation (B 10) is positive if either § > 1 or 

0 < £ < 1 arid p < /? max , where 

Pmax = -\——p2 ‘ (Bll) 

So, for § > 1 flow will accelerate monotonically, while for 0<? < 

1 the momentum will oscillate in the range [0, p max ]. 

According to equation (B10) dp/ds does not depend explicitly 

on the distance s and so its absolute value is the same for the same 
value of p\ therefore for 0<? < 1 the function p(s) is symmetric 
around its maxima and minima and is periodic. 


B2 Ultra-relativistic flow 


For p 1 we can neglect terms of the order 0(1 /p) and higher. 
Then equation (B10) takes the form 


dp\ : 
ds J 


= 2[l + />($-l)] ; 


(B12) 


the solution of equation (B12) is (cf. equation B3 in Fawley et al. 
1977) 


= ~Jls + s 2 . 


(B13) 


If? > 1 the flow is continuously accelerating and its momentum 
at large s will grow as 


?-l 2 

P = — * • 


(B14) 


and the Debye length (also the skin depth, since the characteristic 
velocity is c) of such a plasma is 


c 

J-D.GJ = 

Wp.GJ 


(B6) 


Introducing normalized distance s = x/Xo, gj equation (B4) be- 
comes 


d 2 y = „ Y 
ds 2 v/y 2 - 1 


(B7) 


For ? < 1 the flow is oscillatory and its momentum periodically 
reaches p max (see equation Bll). If 1 — £ 1 and p max 1 

the flow is relativistic almost everywhere. The wavelength of such 
spatial oscillations is twice the distance between the points where 
p = 0 and p reaches its maximum value. From equation (B13) the 
distance where p reaches it maximum is x/2(l — ?) _1 and so the 
period of spatial oscillation is (Beloborodov 2008) 

.50 = 272(1 — tr 1 . 


(B15) 
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B3 Non-relativistic flow 


Stationary space charge limited flow is non-relativistic at the be- 
ginning, close to the NS surface; it can remain non-relativistic if 
f <^C 1 (see equation B 1 1). 

For p 1 we can neglect terms of the order 0(p 3 ) and higher 
and equation (BIO) becomes the Cycloid equation 


f d pV = 2g - p 
\dy J p 


(B16) 


The solution of this equation in parametric form is (Beloborodov 
2008) 


P = ?(1 - COS Wp.Gjt) 


(B17) 


5 = £(<Wp,Gjf - sino>p iG jf) , (B18) 

where the time t is measured in seconds. The period of spatial 
oscillations is then 

s 0 = 2n^ . (B19) 

For small £ the flow is always non-relativistic and equations (B17) 
and (B18) describe it accurately everywhere. For large f equations 
(B17) and (B18) are good approximation near the starting point of 
the flow for all £ and for f < 1 also near periodically repeating 
stagnation points, with small values of p. 

For small s, near the flow’s starting point, the time t can be 
eliminated from (B17) and (B18) and we get (cf. equation 13 in 
Fawley et al. 1977) 

P - ' ? 17 V /3 - (B20) 


APPENDIX C: BOUNDARY CONDITIONS 
FOR MODELLING OF SPACE CHARGE 
LIMITED FLOW 

For modelling of the space charge limited flow in the polar cap one 
needs to make an adequate numerical model for an infinitely large 
pool of particles available at the NS surface in order to correctly 
simulate the space charge limitation condition. This is not a trivial 
task. We found the following procedure works well. 

The calculation domain of the length L is divided in M x equal 
numerical cells (a typical value of M, in our simulations is ~ a few 
thousand). The electric field E and current j are set at cell bound- 
aries i = 0, , M x . For the calculation of the current density we 

use a ID version of the charge conservative algorithm proposed by 
Villasenor & Buneman (1992), when charged particles are repre- 
sented by uniformly charged sheets with the width equal to the cell 
size Ax and the position of a particle is the position of the sheet’s 
centre. The fraction of the sheet passed through the cell boundary 
i during a time-step determines the contribution of the particle into 
the current /,■ at that point. 

At each end of the calculation domain we have one ‘ghost’ cell. 
The outside boundaries of the ghost shells are ‘ghost’ points with 
indexes i = — 1 and i = M x + 1. Equation (12) for the electric 
field is solved for points 0 . . . M x . Particles can move into the ghost 
cells but when their positions are outside of the domain [ — Ajc/2; 
L + Ax/2] and they are moving outwards they do not contribute to 
the current density in the domain anymore and such particles are 
deleted at the end of each time-step. 


x : 


ghost cell 

computational domain 




numerical particle 



Ax 



Figure Cl. Numerical implementation of boundary conditions for space 
charge limited flow. Injection of one numerical particle is shown; particle’s 
centre is marked by a cross. See text for explanations. 


The electric field at the ghost points is set to zero. We solve 
a ID initial value problem - the electric field in the domain is 
calculated from the values of the electric field at the same points 
at the previous time-step and is not coupled to the electric field 
at the ‘ghost’ points. The electric field inside the ghost cells, at 
the particles’ position, is obtained by quadratic interpolation using 
values E-i, Eo, E\ (and Em x -i, Em x , Em x + i) and, therefore, setting 
E at the ghost point to zero (£_ i = E Mx+ 1 = 0) smoothly reduces 
the electric field inside ‘ghost’ cells towards their outer ends. Setting 
the electric field at ghost points at each time-step to the values 
obtained as extrapolation of electric field values near the domain 
boundaries (e.g. using quadratic extrapolation from points 0, 1, 2 
and M x — 2, M x — 1, M x ) or to some non-zero values resulted only 
in higher numerical noise and did not change the system behaviour. 

At the beginning of each time-step we inject certain amount of 
electrons V; n j and equal number of heavy positive charged particles 
(‘ions’) in the first ‘ghost’ cell of our computation domain at the 
position slightly outside the centre of the fist cell x mj = —(Ax/2 
+ 5) (see Fig. Cl). The momentum of each injected particle is 
sampled from a uniform distribution in the interval [— pi„j, p; n j]. In 
this way, we can model finite temperature of the particles on one 
hand and populate the domain with particles more uniformly on 
the other hand, as each injected particle after the first time-step will 
have a slightly different position. If during this time-step the electric 
field inside the ‘ghost’ cell cannot move particles into the interval 
[—Ax/2; L + Ax/2], they do not contribute to the current density 
and will be deleted at the end of the time-step. Depending on the 
value of the electric field either positive or negative particles will 
be ‘sucked’ into the domain. 

We experimented with different values for and found that usu- 

ally after Vi n j exceeds the critical value VV = 2(j m / Qc)(cAt / Ax ) 
by ~ 20-30 per cent computational results stop depending on /V; n j ; 
further increase in /Vmj results in higher numerical noise. 1VV is 
twice the number density of particles with relativistic velocities 
which must be injected at every time-step in order to provide the 
required current density j m \ Q is the charge of a numerical particle 
and c is the speed of light. The factor of 2 accounts for injected 
particles having negative initial momentum - most of them do not 
reach computational domain and are deleted. Having some injected 
particles with negative momenta results in slightly lower numerical 
noise as well as more realistically represents the finite temperature 
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of the NS atmosphere (the equivalent of the warm cathode in the 
analogous vacuum tube and high current beam technologies). The 
computational overhead caused by such particles is negligible, as 
lVj n j is orders of magnitude less than the total number of particles. 
We also experimented with different values for the time-step and 
found that values of At such that Ax/cAt ~ 5 results in relatively 
low level of numerical noise due to discrete events of particle injec- 


tion; smaller At leads to larger numerical overhead, as the stability 
of the leapfrog scheme requires only that At < 0.5 A x/c. 
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